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LSPRAY-V: A Lagrangian Spray Module 


M.S. Raju 

Vantage Partners, LLC 
Brook Park, Ohio 44142 


ABSTRACT 

LSPRAY-V is a Lagrangian spray solver de- 
veloped for application with unstructured grids and 
massively parallel computers. It is mainly designed 
to predict the flow, thermal and transport properties 
of a rapidly vaporizing spray encountered over a wide 
range of operating conditions in modern aircraft en- 
gine development. It could easily be coupled with any 
existing gas-phase flow and/or Monte Carlo Proba- 
bility Density Function (PDF) solvers. The manual 
provides the user with an understanding of various 
models involved in the spray formulation, its code 
structure and solution algorithm, and various other 
issues related to parallelization and its coupling with 
other solvers. With the development of LSPRAY-V, 
we have advanced the state-of-the-art in spray com- 
putations in several important ways. 

• It facilitates the use of both structured & un- 
structured grids and parallel computing and, 
thereby, facilitates large-scale combustor com- 
putations involving complex geometrical config- 
urations. The solver accommodates the use of 
an unstructured mesh with mixed elements of 
either triangular, quadrilateral, and/or tetrahe- 
dral type. 

• In order to deal with modern gas-turbine fuels 
that are mixtures of many compounds, it takes 
into account the modeling of multicomponent 
liquid fuels with variable properties. 

• Various well-established vaporization and atom- 
ization models are incorporated into our spray 
code to cover a wide range of engine operating 
conditions: low to high pressures including su- 
percritical conditions, and superheat conditions 
associated with flash vaporization. The initial 
droplet conditions could be prescribed based on 


either a single-point or multi-point droplet injec- 
tion. The multi-point injection could be in the 
form of a line or circular point injection. The 
initial droplet conditions could be specified en- 
tirely in the form of a table based on known ex- 
perimental data, or some of the initial conditions 
could be calculated based on several widely-used 
droplet-size distribution functions or primary at- 
omization models incorporated into our spray 
code. For low pressures, we have incorporated 
the following primary atomization models: (a) 
sheet breakup, (b) air blast, (c) blob jet, or (d) 
BLS (Boundary-Layer Stripping). We have also 
incorporated a flash induced atomization model 
to cover superheat conditions. The secondary 
droplet breakup following primary atomization 
could be modeled by one of the following mod- 
els: (a) Rayleigh- Taylor, (b) TAB (Taylor Anal- 
ogy Breakup), or (c) ETAB (Enhanced Taylor 
Analogy Breakup). 

• The spray module has a multi-liquid and multi- 
injector capability. 

• Our spray code supports all of the boundary con- 
ditions of the national combustion code includ- 
ing particle movement through very complex pe- 
riodic boundary conditions. Upon impact with 
the wall, a droplet may shatter, rebound, or stick 
to the wall depending on the level and condi- 
tions of collision impact. We have implemented 
several droplet-wall interaction models into our 
spray code to cover a wide range of conditions. 

• Our spray code can be used in the calculation of 
both steady as well as unsteady computations. 

• The spray module is designed in such a way so 
that it could easily be coupled with other CFD 
codes. 
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With the aim of improving the overall solution 
procedure of the national combustion code involving 
sprays, we have made several other relevant contri- 
butions to the gas-side of computations: 

• In order to demonstrate the importance of chem- 
istry/turbulence interactions in the modeling 
of reacting sprays, we have extended the joint 
scalar Monte Carlo PDF (Probability Density 
Function) approach to the modeling of spray cal- 
culations with unstructured grids, and parallel 
computing. 

• In order to account for non-ideal gas behavior as- 
sociated with high-pressure conditions, we have 
completed the implementation of a high-pressure 
EOS (equation of state). Also, we have imple- 
mented several high-pressure corrections that go 
into the calculation of gas-phase transport prop- 
erties. 

The spray solution procedure provided favor- 
able results when applied to the modeling of vari- 
ous reacting/non-reacting flows encountered in gas- 
turbine combustors, stratified-charge rotary combus- 
tion (Wankel) engines, supersonic and pulse detona- 
tion combustion devices. Its use has been demon- 
strated in various NASA projects: the NASA’s fun- 
damental aeronautics program initiative on emissions 
through the Supersonics (SUP) and Subsonic Fixed 
Wing (SFW) project offices, Ultra-Efficient Engine 
Technology (UEET), Pulse Detonation Combustion 
Technology (PDCT), & Rotary Combustion Engine 
Technology Enablement Project (RCETEP). 


N AS A/CR— 20 15-218918 


2 


NOMENCLATURE 


a 

a 0 


A 

B k 

B t 

B 0 

Cfreq 

c 

c D 

c p 

D 

Dab 

d 

d b 

dL 

do 

dparent 

dproduct 

dt 

dO 

Dparent 

Dproduct 

h 

Iqi h 
k 


K b 

K sb 

ki,k 2 

Km 

Ko-K, 


l 

lk,eff 

^kn 

L 


liquid jet radius, m, or 
parent drop radius, m 
initial droplet radius, m 
outward area normal vector 
of the nth surface, m 2 
a constant 

Spalding mass transfer number 
Spalding heat transfer number 
a constant 
a constant 
a constant 
drag coefficient 
specific heat, J/(kg K) 
turbulent diffusion coefficient, m 2 /s 
diffusion coefficient, cm 2 /s 
droplet diameter, m 

droplet diameter after primary breakup, m 

ligament diameter, m 

orifice diameter, m 

parent drop diameter, m 

product drop diameter, m 

time increment, s 

half-cone angle, deg. 

energy contained in the parent 

drop, (kg m 2 )/s 2 

energy contained in the product 

drop, (kg m 2 )/s 2 

specific enthalpy, J/kg, 

liquid sheet thickness, m 

modified Bessel functions of the first kind 

turbulence kinetic energy, m 2 /s 2 , 

wave number, 1/m, or 

thermal conductivity, J/(ms K) 

binary interaction parameter 

a breakup constant, 1/s 

wave number associated with 

sheet breakup, 1/m 

constants 

wave number associated with 
ligament breakup, 1/m 
modified Bessel functions of the 
second kind 

modified wave number, 1/m 

mixture latent heat of evaporation, J /kg 

effective latent heat of evaporation, 

J/kg (defined in Eq. (35)) 
heat of vaporization at normal 
boiling point, J/kg 
liquid sheet breakup length, m 


m 

liquid mass flow rate, kg/s 

rh k 

droplet vaporization rate, kg/s 

rilk, flash 

droplet vaporization rate under 
flash evaporating conditions, kg/s 

rilko 

initial mass flow rate associated 
with kth droplet group, kg/s 

rn k ,t 

droplet vaporization rate due 
to heat transfer, kg/s 

TO 0 

mass of the parent drop, kg 

fh(t) 

mean mass of the product drop distribution, kg 

M 

molecular weight, kg/kg-mole 

M a 

molecular weight of gas excluding 
fuel species, kg/kg-mole 

M f 

molecular weight of fuel, kg/kg-mole 

Mi 

molecular weight of ith 
species, kg/kg-mole 

M s hed 

shed drop mass, kg 

MMD 

mass mean diameter, m 

n k 

number of droplets in kth group 

Tl parent 

drop number in the parent group 

product 

drop number in the product group 

n(t) 

non-dimensional number ( =mo/fh(t )) 

N 

drop number 

N 0 

parent drop number 

N f 

number of surfaces contained in 
a given computational cell 

N p 

total number of computational cells 

Nu 

Nusselt number 

P 

pressure, N/m 2 

Pc 

critical pressure, N/m 2 

Pr 

Prandtl number 

P r 

reduced pressure, P/P c 

P sat 

saturation pressure, N/m 2 

Q 

density ratio {=p g / pi) 

f 

radius (=y/a 3 /r S MR), m 

r k 

droplet radius, m 

r ko 

initial droplet radial location, m 

rsMR 

Sauter mean radius, m 

R u 

gas constant, J/(kg K) 

Re 

Reynolds number 

RND 

random number 

Sh 

Sherwood number 

s k 

droplet radius-squared ( = r 2 ), m 2 

Smlc 

source term contribution from liquid 
exchange to mass conservation, kg/s 

Smle 

source term contribution from liquid 
exchange to energy conservation, J/s 

Smlm 

source term contribution from liquid 
exchange to momentum conservation, kg m/s 2 
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Smis source term contribution from liquid 
exchange to species conservation, kg/s 

SMD Sauter mean diameter, m 
t time, s 

td non-dimensional time (= 2 ((/° ) 

T temperature, K, or 

non-dimensional number (=ZW e°' 5 ) 

Tf, boiling temperature, K 

T n b normal boiling temperature, K 

T c critical temperature, K 

Tfc kth droplet temperature, K 

T r reduced temperature, T/T c 

iid drop deceleration, m/s 2 

Ui ith velocity component, m/s 

Uik ith velocity component of 

kth droplet group, m/s 

U relative velocity between gas and liquid, 
m/s, or relative drop velocity, m/s 
vt product drop velocity component, m/s 

V initial liquid velocity, m/s 

V c volume of the computational cell, m 3 , or 

critical molar volume, m 3 /kg-mole 
V n molar volume at normal pressure, 

m 3 /kg-mole 

We Weber number ( = £aML) 

ibj gas-phase chemical reaction rate, 1/s 

x Cartesian coordinate, m, or drop 

deformation distance, m 

Xi Cartesian coordinate in the ith direction, m 

y Cartesian coordinate, m, or non-dimensional 

deformation parameter (=yf ) 
yj mass fraction of jth species 

x spatial vector 

z Cartesian coordinate, m 

Z compressibility factor (=^p), or 

non-dimensional number (= W ^ l ) 

Z c critical compressibility factor (=AAk) 

a represents a coordinate related to 

a Hill’s Vortex, spray cone rotation 
angle, deg., or a constant 

a s overall heat transfer coefficient, kJ/sm 2o K 

aq a constant 

«2 a constant 

/3 spray cone rotation angle, deg. 

X mole fraction 

5 Dirac-delta function, or 

initial liquid sheet thickness, m 

A p pressure drop in the injector, N/m 2 

A tf local time step used in the flow solver, s 


A t g i global time step used in the spray solver, s 

Atu fuel injection time step, s 

At mi allowable time step in the spray solver, s 
AV computational cell volume, m 3 

e rate of turbulence dissipation, m 2 /s 3 

ej fractional mass evaporating rate of species 

at the droplet surface 
rj wave amplitude, m 

r 0 turbulent diffusion coefficient, kg/ms, or 

a factor (=210 [ Tc pl j^ 6 ) 

A thermal conductivity, J/(ms K), or 
wavelength, m 

A wavelength associated with fl, or 

wavelength associated with 
Rayleigh- Taylor breakup, m 
/i dynamic viscosity, kg/ms 

v kinematic viscosity, m 2 /s 

ui turbulence frequency, 1/s, 

frequency (=^s - ^r), 1/s, 
complex growth rate, 1/s, or 
acentric factor 

fl maximum growth rate, 1/s 

p density, kg/m 3 

pin liquid density (at 1 bar, 273.15 K), kg/m 3 
p r reduced density (=p/ p c ) 

a potential length constant, Angstrom (=0.1 

surface tension, kg/s 2 , or characteristic 
diameter of a molecule in Table 1 
t stress tensor term, kg/ms 2 , or 

characteristic breakup time (=1/A'b), s 
6 void fraction, or 

spray cone angle, deg. 

Subscripts 

A A-th species component 

b breakup conditions 

B B-th species component 
c critical conditions 

/ fuel 

g gas or global 

i i-th coordinate, i-th species 

component, summation index, or 
imaginary component 
inj injector 

j j-th coordinate, j-th species 

component, or summation index 
k k-th droplet group, k-th coordinate, 

summation index, or liquid conditions 
associated with k-th droplet group 
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1 liquid 

L liquid ligament 

m multi-component mixture 

n normal, or nth-face of the 

computational cell 
o initial conditions, 

orifice exit conditions, 
injector initial conditions, or 
oxidizer 

p particle, drop, or conditions 

associated with a grid cell 
r reduced, radial coordinate, real 

component, or a reference value 
RT Rayleigh- Taylor 

s droplet surface, or adjacent 

computational cell 
t time 

x axial or x-coordinate 

y y-coordinate 

2 z-coordinate 

, partial differentiation with respect 

to the variable followed by it 

Superscripts 

mean, or average 

first order differentiation, or 

flow rate 

second order differentiation 
// fluctuations 
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1 INTRODUCTION 

There are many occurrences of sprays in a vari- 
ety of industrial and power applications and materi- 
als processing [1]. A liquid spray is a two phase flow 
with the gas as the continuous phase and the liquid 
as the dispersed phase in the form of droplets or lig- 
aments [1]. The coupling between the two phases oc- 
curs through the exchanges of mass, momentum, and 
energy involving a wide range of thermal, mass, and 
fluid dynamic factors. A number of finite-difference 
formulations have been advanced over the years for 
predicting the flow (mass and momentum) and ther- 
mal properties of a rapidly vaporizing spray. Some 
of the pros and cons of various formulations can be 
found in [1-5]. 

The state of the art in multi-dimensional com- 
bustor modeling, as evidenced by the level of sophis- 
tication employed in terms of modeling and numeri- 
cal accuracy considerations, is also dictated by the 
available computer memory and turnaround times 
afforded by present-day computers. With the aim 
of advancing the current multi-dimensional com- 
putational tools used in the design of advanced 
technology combustors, a solution procedure is de- 
veloped that combines the novelty of the coupled 
CFD/spray/scalar Monte Carlo PDF (Probability 
Density Function) computations on unstructured 
grids with the ability to run on parallel architec- 
tures. In this approach, the mean gas-phase velocity 
and turbulence fields are determined from the solu- 
tion of a conventional CFD method, the scalar fields 
of species and enthalpy from a modeled PDF trans- 
port equation using a Monte Carlo method, and a 
Lagrangian-based dilute spray model is used for the 
liquid-phase representation. 

The gas-turbine combustor flows are often char- 
acterized by a complex interaction between various 
rate-controlling processes associated with turbulent 
transport, mixing, chemical kinetics, evaporation and 
spreading rates of spray, convective and radiative 
heat transfer, among others [10]. The phenomena 
to be modeled as controlled by these processes often 
interact with each other at various disparate time and 
length scales. In particular, turbulence plays an im- 
portant role in determining the rates of mass and heat 
transfer, chemical reactions, and liquid-phase evapo- 
ration in many practical combustion devices. The in- 
fluence of turbulence in a diffusion flame manifests it- 
self in several forms, ranging from the so-called wrin- 
kled or stretched flamelets regime to the distributed 


combustion regime, depending upon how turbulence 
interacts with various flame scales [69-70]. 

Most of the turbulence closure models for reac- 
tive flows have difficulty in treating nonlinear reaction 
rates [69-70]. The use of assumed shape PDF meth- 
ods was found to provide reasonable predictions of 
pattern factors and NO.y emissions at the combustor 
exit [71]. However, their extension to multi-scalar 
chemistry becomes quite intractable. The solution 
procedure based on the modeled joint-composition 
PDF transport equation has an advantage in that it 
treats the nonlinear reaction rates without any ap- 
proximation. This approach holds the promise of 
modeling various important combustion phenomena 
relevant to practical combustion devices such as flame 
extinction and blow-off limits, and unburnt hydrocar- 
bons (UHC), CO, and NOx predictions [71]. 

With the aim of demonstrating the viability of 
the PDF approach to the modeling of practical com- 
bustion flows, we have undertaken the task of extend- 
ing this technique to the modeling of sprays as a part 
of the NCC (National Combustion Code) develop- 
ment program [7-15]. In order to facilitate large-scale 
combustor computations, we have extended our pre- 
vious work on the combined CFD/spray/PDF com- 
putations to parallel computing [10-15]. The use 
of parallel computing offers enormous computational 
power and memory as it can make use of hundreds 
of processors in concert to solve a complex problem. 
The trend towards parallel computing is driven by 
two major developments: the widespread use of dis- 
tributed computing and the recent advancements in 
MPPs (Massively Parallel Processors). The solver is 
designed to be massively parallel and automatically 
scales with the number of available processors. A cur- 
rent status of the use of the parallel computing in tur- 
bulent reacting flows involving sprays, scalar Monte 
Carlo PDF, and unstructured grids was described in 
Ref. [II]. It outlined several numerical techniques 
developed for overcoming some of the high computer 
CPU-time and memory-storage requirements associ- 
ated with the use of Monte Carlo solution methods. 
The parallel performance of both the PDF and CFD 
modules was found to be excellent but the results 
were mixed for the spray computations showing rea- 
sonable performance on massively parallel computers 
like Cray T3D; but its performance was poor on the 
workstation clusters [10]. In order to improve the 
parallel performance of the spray module, two differ- 
ent domain decomposition strategies were developed 
and the results from both strategies were summarized 
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[10-14], 

It is also well known that considerable effort 
usually goes into gridding up of complex gas-turbine 
combustor geometries. In order to allow repre- 
sentation of complex geometries with relative ease, 
we have extended our previous work on the com- 
bined CFD/spray/PDF computations to unstruc- 
tured meshes [11-15], The grid generation time asso- 
ciated with gridding up practical combustor geome- 
tries, which tend to be very complex in shape and 
configuration, could be reduced considerably by mak- 
ing use of existing unstructured grid generators. The 
solver accommodates the use of an unstructured mesh 
with mixed elements: triangular and/or quadrilateral 
for 2D (two-dimensional) geometries and tetrahedral 
for 3D. 

With the development of two computer mod- 
ules, LSPRAY - a Lagrangian spray solver [14] and 
EUPDF - an Eulerian Monte Carlo PDF solver [15], 
we were able to demonstrate the use of the joint 
scalar Monte Carlo PDF method in the modeling of 
complex multidimensional reacting flows (e.g., gas- 
turbine combustor flows) [10-15]. In this manual, we 
only concentrate on providing the details of our spray 
module. However, further details on the application 
of the joint scalar Monte Carlo PDF method to two- 
phase flows could be found elsewhere in [10-13,15]. 

In this manual, we summarize some important 
aspects of our spray formulation without making any 
attempt to provide an in-depth review on the fluid 
dynamic and transport behavior of reacting sprays 
[6-15]. Depending on the nature of the spray, an 
appropriate selection could be made from the choice 
of various well-known spray formulations (multicon- 
tinua, discrete-particle, or probabilistic) based on ei- 
ther a Lagrangian or an Eulerian representation for 
the liquid-phase equations by making use of appro- 
priate droplet sub-grid models. The present solu- 
tion procedure could be used within the context of 
both multicontinua and probabilistic spray formula- 
tions, as it allows for resolution on a scale greater 
than the average spacing between two neighboring 
droplets [1]. For NCC, the adopted choice for the 
gas phase was an Eulerian scheme. The liquid-phase 
equations form a system of hyperbolic equations and 
they could be solved by means of either an Eulerian or 
a Lagrangian representation. A Lagrangian scheme is 
chosen as it reduces the errors associated with numer- 
ical diffusion. The liquid-phase formulation is based 
on various well-established models for droplet drag; 
the vaporization models of a polydisperse spray take 


into account the transient effects associated with the 
droplet internal heating and the forced convection ef- 
fects associated with droplet internal circulation; and 
it employs models for gas-film valid over a wide range 
of low to intermediate droplet Reynolds numbers [7]. 
Our current formulation is applicable for flows with a 
dilute spray approximation where the droplet loading 
is low. The numerical method could be used within 
the context of both steady and unsteady calculations 
[8-13]. Not considered in the present release of the 
code are the effects associated with droplet/shock in- 
teraction and dense spray effects. 

Currently, most of the finite-difference tech- 
niques used for predicting the two-phase flows make 
use of the physics derived from single-component liq- 
uid droplet studies with constant properties. How- 
ever, it is well known that most of the gas-turbine fu- 
els are multicomponent mixtures of many compounds 
with a wide distillation curve [7,18]. The multicom- 
ponent nature of the liquid sprays is becoming evi- 
dent with the increasing need to use jet fuels derived 
from heavier petroleum compounds. The gasifica- 
tion behavior of a multicomponent fuel droplet may 
differ significantly over that of a pure single compo- 
nent fuel droplet [18]. Also, the calculation of the 
variable thermo-transport properties of the liquid- 
mixtures becomes more important at high pressures. 
The flame ignition characteristics (such as the phe- 
nomena associated with flame blow-off and extinction 
conditions) could also be influenced by the nonuni- 
form concentration of the fuels with different volatili- 
ties. However, the importance of the multicomponent 
liquid fuels with variable properties received little at- 
tention in the modeling of comprehensive gas-turbine 
combustor spray computations. With this in mind, 
we have extended the spray formulation to multicom- 
ponent liquid sprays in order to deal with the gas 
turbine fuels that are mixtures of many compounds. 
This implementation also takes into account the ef- 
fect of variable liquid properties. 

In order to reduce some uncertainty associated 
with the specification of initial droplet conditions, 
we have undertaken the task of integrating an at- 
omization module into our spray solution procedure. 
The atomization module was developed by CFDRC 
Inc. [47] in collaboration with the university of Wis- 
consin (UW) [45-46]. Our computational experience 
with the assessment of various atomization models 
was summarized in [59-60,78]. A complete descrip- 
tion of all the primary atomization as well as the 
secondary droplet breakup models contained in the 
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CFDRC/UW atomization module is provided in the 
appendix. 

An understanding of droplet behavior under su- 
perheat conditions is identified to be a topic of im- 
portance in the design of modern supersonic engines, 
scramjet, and ramjet afterburners [61-64, & 74]. It is 
because (1) higher heat-load requirements may lead 
to the use of the same fuel as a coolant, (2) under 
some engine operating conditions the nozzles are re- 
quired to operate under low backpressures, and (3) 
under normal gas-turbine operating conditions, it is 
estimated that a small fraction of the liquid fuel 
may be released by flash boiling [62]. There are 
some reported incidences of engine damage due to 
flash atomization in gasoline direct-injection internal 
combustion engines [62]. Because of its importance 
in some NASA-supported future engine development 
projects, we have undertaken an assessment study to 
establish the baseline accuracy of various existing va- 
porization and atomization models used in the calcu- 
lation of a superheated spray. As a part of this effort, 
we have incorporated a solution procedure based on 
the modeling approach developed by [62-64, & 74], 
Our computational experience with the assessment 
of this modeling approach is summarized in [78-79]. 

There is a need to understand the droplet be- 
havior under supercritical pressures as some mod- 
ern gas-turbine combustors tend to operate under 
increasingly high pressures [19-20, & 81-86]. Based 
on some important aspects of the modeling approach 
derived from [20, & 83-86], we have implemented 
a high-pressure droplet vaporization model into our 
spray code. As a part of this effort, we have also 
implemented the Peng-Robinson equation of state 
following the approach of [16-17, & 28], the high- 
pressure corrections to the gas-phase transport prop- 
erties from [16, 19, & 29-31], and the high-pressure 
corrections to the liquid-phase transport properties 
from [16-18]. 

One important aspect fo spray modeling is how 
to deal with the collision dynamics of a droplet im- 
pact with the wall. Upon its impact with the wall, a 
droplet may shatter, rebound, or stick to the wall. In 
order to deal with several possible outcomes, we have 
incorporated the modelng approach of [80-81] based 
on the coding partially received from CFDRC. 

Some other important aspects of the spray mod- 
ule are summarized below: 

• An efficient particle tracking algorithm was de- 
veloped and implemented into the Lagrangian 


spray solver in order to facilitate particle move- 
ment in an unstructured grid of mixed elements. 

• LSPRAY-V is currently coupled with an un- 
structured flow (CFD) solver of NCC [21-23], 
and an Eulerian-based Monte Carlo probability 
density function solver - EUPDF [15] , which were 
selected to be as the integral components of the 
NCC cluster of modules. EUPDF provides the 
solution for the scalar fields of species and en- 
thalpy from a modeled PDF transport equation 
using a Monte Carlo method. 

• The spray solver receives the mean velocity and 
turbulence fields from the flow solver. The so- 
lution for the scalar (energy and species) fields 
could be provided by means of either a conven- 
tional CFD solver or a Monte Carlo PDF solver 
depending on the choice of the solver. 

• The spray solver supplies the spray source-term 
contributions arising from the exchanges of mass, 
momentum and energy with the liquid-phase 
to the flow solvers (CFD and/or Monte Carlo 
PDF). This information could be used in either 
conservative or non-conservative finite-difference 
formulations of the gas phase equations. 

The furnished code demonstrates the success- 
ful methods used for parallelization and coupling of 
the spray to the flow code. First, complete details of 
the spray solution procedure is presented along with 
several other numerical issues related to the coupling 
between the CFD, LSPRAY-V, and EUPDF solvers. 
It is followed by a brief description of the combined 
parallel performance of the three solvers (CFD, EU- 
PDF, and LSPRAY-V) along with a brief summary 
of the validation cases. 

2 GOVERNING EQUATIONS FOR GAS 
PHASE 

Here, we summarize the governing gas-phase conser- 
vation equations in Eulerian coordinates [1]. This 
is done for the purpose of identifying the interphase 
source terms arising from the exchanges of mass, mo- 
mentum, and energy with the liquid phase. They are 
valid for a dilute spray with a void fraction of the 
gas, 9 1 close to unity. The void fraction is defined as 
the ratio of the equivalent volume of gas to a given 
volume of a gas and liquid mixture. 

The conservation of the mass leads to: 
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\pV c ],t + \pV c Ui\ tXi = SmZc = ^2 Uk ™ k ( 1 ) 

k 

The source term is given as a summation over dif- 
ferent classes of droplets. Each class represents the 
average properties of a different polydisperse group of 
droplets. Here, n k represents the number of droplets 
in a given class and m k represents the corresponding 
mass vaporization rate. 

For the conservation of the jth species, we have: 

\pV c yj],t + IpVcV+yjltXi ~ [pV c Dyj iX .] iXi - pV c Wj = 

S-mls = ^ ^ Tlk til 'k ( 2 ) 

k 

where 


Similarly, the energy conservation has a source term 
contribution from the added liquid fuel vapor and an 
additional contribution associated with the effective 
latent heat of vaporization which accounts for the 
heat flux (loss or gain) to the droplet interior from 
the ambient. 

The main purpose of the spray solver is to cal- 
culate the source terms arising from the exchanges 
of the mass, momentum, and energy. In the case of 
NCC, it provides the calculated source terms to both 
the CFD and Monte Carlo PDF solvers. 

3 HIGH PRESSURE EQUATION OF 
STATE 

In order to calculate the high pressure gas be- 
havior, the Peng-Robinson EOS (Equation-Of-State) 
is employed for a multi-component mixture in the fol- 
lowing form [16-17,28]: 



0 and 6j 

3 


1 


P = 


RT 


V-b 


m 


dm 

U 2 + 2 b m V - bl 


For the species conservation, the source term con- where 
tains an additional variable, e,-, which is defined as n 

the fractional vaporization rate for species j. 

For the momentum conservation, we have: h 

u m 


E i Ej yiyj(cnaj ) 1/2 ( l - %), 
Ei?/A, 


( 5 ) 


[pV c Ui\, t + [, pV c UiUj \ Xj + [pV c ], Xi - 

\0V c Tij] } Xj [(1 ) V c Ti-ij\ x .j — S m i m = 


bi = 


0.07780 RT ic 


CLi — 


0.45724 R Z T- 

Pic 


l [l + fiU 1 ~ T ir /2 )] 2 , 


y. n k m k u ki ~y ^ p k r k n k u kijt (3) 

k k 

where the shear stress r,;j in Eq. (3) is given by: 

7 'ij — p\ni, X j T Uj^xi\ T^pdijUi^Xj 

For the momentum conservation, the first source 
term represents the momentum associated with liquid 
fuel vapor and the second represents the momentum 
change associated with droplet drag. 

For the energy conservation, we have: 


Ti r = T /Ti c , 

f iu = 0.37464 + 1.54226a;* - 0.26992a;* 2 , 

Wj is known as the acentric factor of the molecules 
which is a measure of non-sphericity of the molecules, 
and kij is known as the binary interaction coefficient. 
However, the Peng-Robinson EOS is rewritten as a 
cubic EOS in terms of the compressibility factor, Z 
(= ^), before it is solved: 

Z 3 - (1 - B*)Z 2 + (A* - 2 B* - 3 B* 2 )Z 


= Smle ~ ^ ' kl k lil k ^h s (4) 


—A*B* + B* 2 + B* 3 = 0 


( 6 ) 


where 


A* = an d B* = 
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Table 1. Physical constants. 

Species 

T n b 

T c 

Pc 

Pc 

l kn 

v c 

v n 

(jJ 

a 

a / k 


(K) 

(K) 

(atm) 

(Kg/m 3 ) 

(KJ/kg) 

(cm 3 /g-mole) 

(cm 3 /g-mole) 


(A °) 

(K) 

CqHh 

341.9 

507.4 

30.0 

660.0 

334.8 

370.0 

140.06 

0.296 

5.949 

399.3 

C-jHiq 

371.6 

540.2 

27.0 

682.0 

316.3 

432.0 

162.00 

0.351 

6.297 

419.031 

CgHis 

398.82 

568.8 

24.6 

718.5 

301.3 

492.0 

188.8 

0.394 

6.62 

488.15 

C10H22 

447.3 

617.6 

20.8 

728.3 

276.1 

603.0 

233.68 

0.490 

7.16 

540.06 

Cl 2^26 

489.5 

658.3 

18.0 

748.0 

256.3 

713.0 

278.54 

0.562 

7.655 

583.68 

C14H30 

526.7 

694.0 

16.0 

763.0 

240.1 

830.0 

326.62 

0.679 

8.067 

629.08 

n 2 

77.4 

126.2 

33.9 

807.1 

197.6 

90.1 

31.87 

0.039 

3.681 

91.5 

0 2 

90.2 

154.6 

50.4 

1135.7 

212.3 

74.4 

26.08 

0.025 

3.433 

113.0 

CO2 

00.0 

304.1 

73.8 

000.0 

000.0 

94.0 

33.32 

0.239 

3.996 

190.0 

H 2 0 

373.2 

647.3 

221.2 

958.1 

2257.2 

57.1 

19.76 

0.344 

2.641 

809.1 


Table 2. Binary Interaction Parameters, kij. 


C7H1Q 

(j=l) 

0 2 

0=2) 

n 2 

0=3) 

C0 2 

0=4) 

h 2 o 

0=5) 

CV - Hl 6 

(i=l) 

0.0000 

0.1321 

0.1440 

0.1000 

0.1484 

0 2 

(i=2) 

0.1321 

0.0000 

-0.0119 

-0.0289 

0.0910 

^ > 
II ^ 

0.1440 

-0.0119 

0.0000 

-0.0170 

0.1030 

C0 2 

(i=4) 

0.1000 

-0.0289 

-0.0170 

0.0000 

0.1200 

h 2 o 

(i=5) 

0.1484 

0.0910 

0.1030 

0.1200 

0.0000 
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We have chosen the Peng-Robinson EOS be- 
cause of its simplicity and, more importantly, it 
proved to be very useful in the supercritical droplet 
vaporization studies of [19-20]. 


(A - A n )r Z\ = 1.2210~ 2 [ea;p(0.535p r .) - 1] (8) 

when p r < 0.5 


Table 1 lists various physical constants for some 
of the species that we found to be of interest in our 
spray computations. It contains values for the boiling 
temperature at normal pressure, the critical tempera- 
ture, the critical pressure, the critical density, the la- 
tent heat of vaporization at normal pressure, the crit- 
ical molar volume, the molar volume at normal pres- 
sure, and the acentric factor of the molecules. Most 
of this data is collected from [16-18] and is useful in 
both the evaluation of the Peng-Robinson EOS and 
the calculation of various other variable properties. 
And Table 2 lists the binary interaction parameters, 
kij , used in a calculation for the multi-component 
mixture of n-heptane, 0 2 , AT 2 , C 0 2 , & H 2 0. While 
the data for most of the binary pairs was obtained 
from various reference books, some of the missing 
data, however, is replaced with known data found for 
other binary pairs of molecules with similar molecular 
weights. 

4 HIGH PRESSURE CORRECTIONS FOR 
GAS-PHASE TRANSPORT PROPERTIES 

The effect of pressure on the viscosity of pure 
gases is determined by making use of the Reichenberg 
method [16,19,29]: 


(A - \n)TZ% = 1.1410 _2 [ea;p(0.67p r ) - 1.069] (9) 
when 0.5 < p r < 2.0 


(A - \ n )TZl = 2.6010~ 2 [exp(1.155p r ) + 2.016] (10) 
when 2.0 < p r < 2.8 


where A is in W/(m.K), T = 210 [ T °pl ] 1 ^ 6 , Z c is the 
critical compressibility, and p r is the reduced density 
p/Pc = Vo/V. 

The effect of pressure & temperature on diffu- 
sion coefficients in gases is determined by making use 
of Takahashi correlation [16,19,31]: 


_D_abP_ 

( D ab P)+ 


f(T r ,P r ) 


( 11 ) 


where Dab = diffusion coefficient, cm 2 /s, P = 
pressure, bar, the superscript + indicates that low- 
pressure values are to be used, and the reduced pres- 
sure and temperature in the above equation are cal- 
culated as follows: 


T r 

T cAB 

Pr 

PcAB 


T 

T c AB 


VaTqA 

p 

PcAB 


VaPcA 


+ VbT c b 
+ VbPcB 


P_ 

Pn 


Qv 


A P 


, 3/2 


B v P r + (1 +C v P r Dv )~ 1 


( 7 ) 


where 


A v = ^exp(a 2 T“) B v = A v (/3iT r - /3 2 ) 
C v = £exp{ l2 Tc) D v = ±ex P (\ 2 T?) 


a\ = 1.982510” 03 
Pi = 1.6552 
7i = 0.1319 
Ai = 2.9496 


a 2 = 5.2683 
P 2 = 1.2760 
y 2 = 3.7035 
A 2 = 2.9190 


a = -0.5767 

c = -79.8678 
d = —16.6169 


and Q v = 1.0 for non-polar molecules. 

The effect of pressure on the thermal conduc- 
tivity of pure gases is determined by making use of 
the Stiel & Thodos method [16,19,30]: 


This correlation is valid for a system involving 
a trace solute diffusing in a supercritical fluid. It is 
shown to provide satisfactory results but it is based 
on a very limited set of available experimental data. 
The graphical representation of the Takahashi func- 
tion is given in Fig. 1. 

It is noteworthy that at low pressures, the poly- 
nomial fits for the variable thermodynamic properties 
are taken from the data set compiled by McBride et cl 
[24]. The transport properties involving the thermal 
conductivity and molecular viscosity for individual 
species is estimated based on the Chapman-Enskog 
collision theory [25-27]. And Wilke’s formulae is used 
to determine the properties of mixture [25-27]. The 
binary-diffusion coefficients are determined based on 
the Chapman-Enskog theory and the Lennard-Jones 
potential [25-27]. 
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5 LIQUID-PHASE EQUATIONS 



0.0 1.0 2.0 3.0 4.0 5.0 6.0 

Reduced pressure, P r 


Figure 1 Takahashi correlation showing the 
variation of the binary diffusion coefficient 
versus reduced pressure at different reduced 
temperatures. 


Here, we summarize the governing equations for 
the liquid-phase based on a Lagrangian formulation 
where the equations for particle position and velocity 
are described by a set of ordinary differential equa- 
tions. For the particle position of the kth droplet 
group, we have: 


dx^ k 

dt 


— y*ik 


For the droplet velocity: 


( 12 ) 


duik 

dt 


3 C D[ x gs Re k 

16 p k r\ 


[U ig +U' g ~ U ik \ 


where 


C D = 


24 

Re k 



(13) 


(14) 


According to Yuen and Chen [32] , the droplet drag for 
a vaporizing droplet could be calculated by a solid- 
sphere drag correlation but suggested using a correc- 
tion in the evaluation of the droplet Reynolds number 
for the average viscosity based on a 1/3 weighting rule 
as given by Eq. (46). The droplet Reynolds number 
is given by: 


Re k = 2 ? -^ [(u g + u' g - u k ) . ( u g + u' g - u k )\ 1 ' 

Pgs 

( 15 ) 


By following the approach taken from KIVA 
[33], a gas turbulence velocity, u' g , is added to the 
local mean gas velocity when calculating the droplet 
drag and vaporization rates. The gas velocity fluctu- 
ations is calculated by randomly sampling a Gaussian 
distribution with mean square deviation, 2/3 k. The 
Gaussian is given by 


G(u' g ) = (4/37rfc) 3/,2 exp[— 3|Mg| 2 /4fc] (16) 

The gas fluctuating component is calculated once at 
every turbulence interaction time, t tur , and is oth- 
erwise held constant [33]. The correlation time is 
taken to be the minimum of either the eddy time or 
the transit time taken by the droplet to traverse the 
eddy. It is given by 


• 

ttur — TTltTL I , Ctt 


fc 3 / 2 


v e 


e \u g + u' g -u k \ 


(17) 


where Ctt is an empirical constant with a value of 
0.16432. 

The liquid mixture density, p k . in Eq. (13) is 
given by [16,18]: 


p k — ) y k i,p k i (18) 

and the individual component liquid density is given 
by: 


M k i 

Pki — -it 
V k i 

where the molar volume, V k i, is 

V ki = V C1 (0.29056 - 0.8775wi) c ’ 

and 


(19) 


(20) 


± Cl J 

The droplet regression rate is determined from 
one of three different correlations depending upon the 
droplet Reynolds number range. The first correlation 
as given by Eq. (21) is based on a gas- film analysis 
developed by Tong & Sirignano [34]. It is based on a 
a combination of stagnation and flat-plate boundary- 
layer analysis and is valid for Reynolds numbers in 
the intermediate range. The last two correlations as 
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given by Eqs. (22) & (23) are taken from Clift et 
al [35] and are valid when Re k < 20. In fact, when 
the droplet Reynolds number goes to zero, Eq. (23) 
becomes identical to the droplet regression rate of a 
vaporizing droplet in a quiescent medium [27]. 


ds k 

dt 



r 2 -I 1 / 2 

-Re h f(Bk) 

7 r 


if Re k > 20 


(21) 


ds k 

dt 


Pa D g 

Pk 


1 + (1 + Refe) 1 ^ 3 


Re° k 077 


ln( 1 


B k ) 

(22) 


if 1 < Re k < 20 


ds k 

dt 


Pg D g 

Pk 


1 + (1 + Refc) 1 ^ 3 ln( 1 + Rfc) 


(23) 


dT k _ A, 

dt Cpipirl 


\ d 2 T k 


da 2 


+ (1 + C7(i)a) 


azy 

9a 


( 27 ) 


where 


C(t) 


3 

17 


(dpi pi 

A; 



(28) 


where a represents the coordinate normal to the 
stream-surface of a Hill’s Vortex in the circulating 
fluid, and C(t) represents a nondimensional form of 
the droplet regression rate. The initial and boundary 
conditions for Eq. (27) are given by: 


a 


a = 1, 


t = 
0, 


m 

da 


t injection ? 

£ 

II 

O 


dT k 

1 

CplPl 

r 2 

dT k 

da 

17 

A; 

r k 

dt 


3 Pk 



ds k 


32 A; 


~ Ik 

dt 


(29) 

(30) 


(31) 


if Re k < 1 

where B k is the Spalding mass transfer number and 
is given by: 


B k 


( Vfs ~ Vf) 
(1 - Vfs) 


(24) 


and yf s is given as a summation over the fuel-species 
mass fractions at the droplet interface: 


vfs = y™ 

i 


(25) 


where a = 0 refers to the vortex center, and a = 1 
refers to the droplet surface, and the mixture latent 
heat of vaporization, l k , is given by 


h = tjhi (32) 

i 

and the individual component latent heat of vapor- 
ization, l k i , is given by [16,18]: 


l 


ki 



T ci ~T k \Q- 38 
T C i — T bi ) 


(33) 


and the droplet boiling temperature is given by 


and yf is given as a summation over the fuel-species 
mass fractions in the ambient: 

Vf = Yl Vfi ( 26 ) 

i 

and the function f(B k ) is similar to that of a Blassius 
function [1, 36] and is obtained from an analysis simi- 
lar to that of Emmon’s boundary-layer flow over a flat 
plate with blowing [36]. The range of validity of this 
function was extended in Raju and Sirignano [7] to 
consider the effects associated with a boundary-layer 
flow with suction. 

The internal droplet temperature is determined 
based on a vortex model [34] . The governing equation 
for the internal droplet temperature is given by: 


Tbi = 


IkinMi/ R u 


IkinMi / ( Rutbni ) - ln(P) 


(34) 


and, finally, the effective latent heat of vaporization, 
l k ,effi is defined as: 


h,eff = h + 47T 


TOfc 


(dn 

\ dr 


S 


(35) 


It is a very useful parameter as it represents the total 
energy loss associated with the latent heat of vapor- 
ization in addition to the the heat loss to the droplet 
interior. l k ,eff is calculated by means of the following 
relationship [18]: 


, _ C p {T g ~T ks ) 

k,eff (i + B k y/ Le - i 


(36) 
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Similar to the internal thermal transport, the 
internal mass transport of a multi-component fuel is 
given by: 


dyki _ „Dk 


a 


d 2 yki 

da 2 


+ (1 + C(t)a) 


dyki 

da 


(37) 


The initial and boundary conditions for Eq. (37) are 
given by: 


t — tinjectiom Vki — 2/fcz,o 


(38) 


a = 0, 


dyki 

da 


1 

17 



dyki 

dt 


(39) 


a = 1, 


dyki 

da 


3_ J_ , _ , dsk_ 

32 D k [Vkis 6ii dt 


(40) 


By knowing the mass fractions of the liquid 
species at the droplet surface, the corresponding mole 
fractions are determined by 


yiks /Mi 

^ ~ Yi Vyks/Mi 


(41) 


At the droplet interface, the mole fractions of 
the gas species are obtained by means of Raoult’s 
law: 


%is — pXiksPis (42) 

where the partial pressure, Pi S , is determined by 
means of the Clausius-Clapeyron relationship: 


P is = exp 


^ki 

( i 

__LY 

R~ u 

KTbi 

Tks ) . 


(43) 


Then, the corresponding mass fractions in the gas- 
phase at the droplet interface are given by 


the temperature and composition as determined by 
the following one-third rule: 

1 2 

1 fiavg = + g 4>ks (46) 

The correlations for the gas-phase thermody- 
namic and transport properties are described in Sec- 
tion 4. In a similar way, the liquid-phase thermody- 
namic and transport properties are determined based 
on the correlations described in the next section. 

6 LOW PRESSURE LIQUID 
THERMODYNAMIC & TRANSPORT 
PROPERTIES 

The specific heat at constant pressure, C p i , 
thermal conductivity, A / , and viscosity, pi , are evalu- 
ated by means of the following expressions: 

Cpi — C-pio + CpiiT + CpnT 2 + Cpi^T 3 + C p i^T a 

(47) 


Ai = A;o + A n T + A l2 T 2 + A l3 T 3 + A /4 P 4 + A l5 T 5 

(48) 

In pi = pio + pn/T + p l2 T + pi 3 T 2 (49) 

where C p i is in J/(kg K), pi in (/^PA s), and A; in 
(W/m K). 

Tables 3-5 provide the polynomial constants 
used in Eqs. (47)-(49) for some of the species that 
we found to be of interest in our spray computations. 
Table 3 provides the constants for the liquid specific 
heat, C p i , Table 4 for the liquid thermal conductivity, 
Ap and Table 5 for the liquid molecular viscosity, /q. 
These tables are compiled with the data taken mostly 
from the references of [16,18]. 

The binary diffusion coefficient, Dij , is evalu- 
ated as follows [16]: 


v . = XisMi (44) 

M aiX ~ Yi X is) + Yi M iXis 

where M a is the molecular weight of the gas excluding 
fuel vapor. And the fractional mass vaporization rate 
of liquid species, , is given by 

i (a \ d is ~ dp (a k'i 

£i = Vis + (1 - Vfs) (45) 

Vfs ~ yf 

It is noteworthy that the thermodynamic and 
transport properties at the gas film are calculated at 


— 


K dif T 


(50) 


where Kdif = 8.210~ s (l + [^y 2 -] 2 ^ 3 )- One should be 
careful in using this approximation as it is based on 
a scarce set of available experimental data. 

The specific heat for a multicomponent mixture 
is given by: 


n 

c P m = Y. y* c P' ( 51 ) 

i= 1 
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Table 3. Polynomial constants for liquid specific heat. 

Fuel 

Cpl 0 

Cpll 

Cpl 2 

Cpl 3 

CplA 

CqHi^ 

2.4169 

-5.9866e-03 

2.0959e-05 

-8.4489e-09 

0.0 

C7-H16 

4.8227 

-3.6980e-02 

1.6777e-04 

-3.0987e-07 

2.2081e-10 

Cs-ffis 

9.2189 

-8.8314e-02 

3.8869e-04 

-7.2539e-07 

5.0776e-10 

C 10 H 22 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

C 12 H 26 

4.7900 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 

C14H30 

4.7991 

-2.8643e-02 

9.3619e-05 

-8.9516e-08 

0.0 



Table 4 

Polynomial constants for liquid thermal conductivity. 


Fuel 

Mo 

Azi 

M 2 

M3 

M4 

Mo 

CqH \4 

0.37078 

-5.4313e-03 

4.628e-05 

-1.8002e-07 

3.2243e-10 

-2.1832e-13 

C 7 H 1 Q 

0.13236 

9.4441e-04 

-6.588d-06 

1.4617e-08 

-1.1244e-ll 

0.0 

Cgili8 

0.25652 

-7.5401e-04 

1.5872e-06 

-1.6795e-09 

-1.3375e-16 

0.0 

C 10 H 22 

0.22179 

-2.3699e-04 

-6.94e-07 

2.0415e-09 

-1.5741e-12 

0.0 

C 12 H 2 S 

0.17609 

4.2463e-05 

-7.4467e-07 

6.9446e-10 

0.0 

0.0 

C 14 H 30 

0.18801 

-9.1399e-05 

-2.1464e-07 

1.1655e-10 

0.0 

0.0 

n 2 

-2.629E-1 

-1.545E-3 

-9.450E-7 

0.0 

0.0 

0.0 

o 2 

2.444E-1 

-8.813E-4 

-2.023E-6 

0.0 

0.0 

0.0 

co 2 

4.070E-1 

-8.438E-4 

-9.626E-7 

0.0 

0.0 

0.0 

H 2 o 

-3.838E-1 

5.254E-3 

-6.369E-6 

0.0 

0.0 

0.0 


Table 5 

Polynomial constants for liquid molecular viscosity. 

Fuel 

[Mo 

Mu 

M12 

Ml 3 

C 5 H 14 

-4.034E+00 

8.354E+02 

0.0000000 

0.0000000 

C 7 H 16 

-4.325E+00 

1.006E+03 

0.0000000 

0.0000000 

Cs-His 

-4.333E+00 

1.091E+03 

0.0000000 

0.0000000 

C 10 H 22 

-4.460E+00 

1.286E+03 

0.0000000 

0.0000000 

C 12 H 26 

-4.562E+00 

1.454E+03 

0.0000000 

0.0000000 

C 14 H 30 

-4.615E+00 

1.588E+03 

0.0000000 

0.0000000 

n 2 

-2.795E+01 

8.660E+02 

2.763E-01 

-1.084E-03 

0 2 

-4.771E+00 

2.146E+02 

1.389E+02 

-6.255E-05 

C0 2 

-3.097E+00 

4.886E+01 

2.381E-02 

-7.840E-05 

h 2 o 

-2.471E+01 

4.209E+03 

4.527E-02 

-3.376E-05 
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8 VAPORIZATION MODELING OF A 
SUPERHEATED DROPLET 


And the thermal conductivity of a multicomponent 
mixture is calculated by means of the Li method [16]. 

n n 

A m = ^ ^ ] c/)i(/>jAij (52) 

i = 1 1=1 

where 

\ij = 2(Aj 1 + A 1 ) _1 

a. - xiYi 


where Xi is the mole fraction of the species i, <j>i is a 
volume fraction of the ith species, and V t is the molar 
volume of the pure fluid. 

7 HIGH-PRESSURE LIQUID TRANSPORT 
PROPERTIES 

The high-pressure correction for thermal con- 
ductivity is calculated by means of a correlation taken 
from Refs. [16-18]. 

Wx = 1.0 + (P r /31.6 + T 2 (1.0 - exp(—0.2P r )))/Q\, P r < 12.0, 

(53) 

Wx = 1.0 + {P° T )( 0.0667 + 0.1597 T?)/Qx,P r > 12.0 (54) 

where Q\ = 8.6514 - 3.78 T r - T r 2 . 

The high pressure correction for viscosity is also 
calculated by means of a correlation taken from Refs. 
[16-18]. 


II 

£ 

(1.0 

+ R A1 (0.47214AP r .) A ' J )/(1.0 + C>AP r ) 

(55) 

where 

A P r 

= 

(P~Psat)/Pc 

A P 

= 

0.999 - 4.442 E - 04/(T r -0 03877 - 0.999) 

Bp 

= 

0.3257/(1.0039 - T 2 ' 673 ) 0 ' 2906 _ 0.2086 

Cp 

= 

— 0.0792 + 2.16167V - 13.404T 2 + 44.1706T 3 

— 84.829 T* + 96.1209T 3 - 59.8127T 3 + 15.6719TJ 


It is noteworthy that in calculating the prop- 
erties of a multi-component mixture that we have 
adopted the same kind of mixing rules as described 
in the previous section. 


Flashing phenomena refers to a process that is 
in thermodynamic non-equilibrium when a liquid is 
superheated [72-73]. The reasons for its occurrence 
are mainly two-fold [72-73]: (1) a liquid fuel can be 
heated to a temperature above its saturation temper- 
ature, and (2) when a liquid is depressurized rapidly 
it can lead to flashing as the thermal inertia tends to 
maintain its bulk internal temperature above the sat- 
uration temperature. Although flash evaporation is 
considered to be detrimental to engine performance 
under normal circumstances, it can have some poten- 
tial benefits. It is known to produce a fine spray with 
enhanced atomization, increase effective spray cone 
angle, and decrease spray penetration [61]. 

An understanding of flash injection is of im- 
portance in some applications involving scramjet and 
ramjet afterburners because the same liquid fuel is of- 
ten used as a coolant. Also, the engine conditions are 
such that nozzles operate at low back pressures with 
supersonic outflow [61]. The objective of our work is 
to establish a baseline accuracy for existing atomiza- 
tion and vaporization models used in the prediction of 
a superheated spray by undertaking a critical review 
of existing experimental data available in the litera- 
ture for validation. This work is funded through the 
supersonics (SUP) and subsonic fixed wing (SFW) 
project office initiatives on high altitude emissions of 
the NASA’s fundamental aeronautics program. 

In what follows we first provide some details of 
the two superheat vaporization models that we have 
implemented into NCC. It is followed by other model- 
ing considerations that need to be taken into account 
as described in Sections 8.3 & 8.4. 

8.1 Superheat Vaporization Model of Zuo, Gomes, & 
Rutland [62], and Schmehl & Steelant [63-64] 

It is based on an extension of the classical D 2 - 
theory. In the classical evaporation model, the ther- 
mal energy needed for evaporation is mostly furnished 
by the external heat transfer from the surrounding 
gas. Under superheat conditions, the characteristic 
vaporization time resulting from the external heat 
transfer from the surrounding gas is of the same order 
of magnitude as that resulting from the flash evap- 
oration. The energy needed for vaporization at the 
droplet surface is partly provided by the superheat 
energy stored within the droplet and it is controlled 
by the droplet internal heat transfer. The modeling 
approach differs from the classical droplet vaporiza- 
tion models in three important ways: (1) under su- 
perheat conditions, the droplet surface mass fraction, 
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Y fs , approaches unity as the droplet surface temper- 
ature is maintained at the corresponding liquid boil- 
ing temperature; (2) under superheat conditions, all 
the external heat transfer from the surrounding gas 
is made available to the vaporization process with no 
apparent increase in the droplet surface temperature; 
and (3) the flow of fuel vapor imparted by flash va- 
porization partly counterbalances the flow generated 
by external heat transfer and may significantly reduce 
the energy transferred from the surrounding gas. 

Based on the governing equations of conserva- 
tion for an isolated spherically symmetric droplet, 
Zuo et al [62] and Schmehl and Steelant [63-64] 
showed that the total evaporation rate, ft i kl can be 
calculated as 


m k = rh k , flash + rh k ,t (56) 

where the flash boiled vaporization rate, frik, flash , is 
given by 


A 2 (T k ~T b ) 

mk, flash = 4? ir k a s (57) 

where T k is the internal droplet temperature and the 
overall heat transfer coefficient, a s (= kJ/s to 2 °K) 
is given by the Adachi correlation [65] : 

= 0.76(T fc - Tfo) 0,26 (0 < T kB -T h < 5) 

a B = 0.027(Tfc - T b ) 233 (5 < T kB - T b < 25) (58) 

= 13.8(T fc - T,,) 0 ' 39 ( T ks -T b > 25) 


The Adachi correlation is valid over a wide range of 
superheat conditions. 

The vaporization rate due to external heat 
transfer, fhk,t , in Eq. (56) is given by 


"ife.t = 271-1-*.- 


Nu 


i j ^ k, flash 


ln[ 1 + (1 + 


ITLk , / lash 


)Bt] (59) 


where the Spalding heat transfer number, B tl and the 
Nusselt number, Nu , are given by 


B t 


C p (T g - T ks ) 

l k,eff 


(60) 


Nu = 2(1 + 0.3i?e 1/2 Pry 3 ) (61) 

The corresponding droplet regression rate, 
is given by 


This model is valid over an entire range of su- 
perheat conditions as long as there is some amount of 
superheat energy available within the droplet (T k > 
T b ). 


8.2 Superheat Vaporization model of Lee et al. [74] 

This modeling approach is also based on an ex- 
tension of the classical P 2 -theory but it differs from 
the previous model of [62, 63-64] in several important 
ways: (1) neglects the effects of internal nucleation 
and convection, (2) instead of invoking the Adachi 
correlation for the flash induced vaporization rate, 
the heat loss from the droplet interior is modeled by 
qi = Ti " teri ° r ~ Ts ( P °°) , (3) the formulation is valid 
only when T, x > Tf,„f,(P 00 ), which means for posi- 
tive Spalding heat transfer numbers. 

With the assumptions employed, Lee et al. [74] 
arrived at the following equation for the evaporation 
rate, m v . 


m v — s = fh v L — qi (63) 

e m v _ 1 

where 

f = C P gT/L r , fhy = (47rr^p s u)/(47rA g r fcit=0 /C'pg), 
L = Lk/L r . The solution for the above equation is 
obtained by employing the standard Newton’s itera- 
tion. 

It is noteworthy that because of the assump- 
tions invoked in this model, the applicability of this 
model may be limited in several important ways: 
(1) our experience shows that the assumption of 
Too > T bub (Poo) becomes too restrictive in most 
non-reactive calculations and also in some regions of 
most reactive spray calculations, and (2) the qi term 
is modeled based on the assumption that the overall 
droplet life is very short (sub millisecond range) but 
in most moderate superheat conditions the droplet 
lifetimes are on par with those vaporizing under non- 
superheat conditions. So its use is likely to be re- 
stricted for the modeling of smaller droplets produced 
after flash induced atomization. 

The solution for Eq. (63) and all the JP8 fluid 
properties needed in this calculation were obtained 
by the computer coding provided by UTRC. 

8.3 Vaporization Model Valid Under Non-Superheat 
Conditions 


dsk _ frik 
dt 2itr k pi 


Under moderate initial superheat conditions, 
only a fraction of the vaporization takes place under 
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9 SUPERCRITICAL DROPLET 
VAPORIZATION 


superheat conditions ( T k > T b ) and the remainder 
takes place under non-superheat evaporating condi- 
tions (Tfc < Tb). So there is a need to revert back to 
a vaporization model valid under stable evaporating 
conditions when the internal droplet temperature ap- 
proaches the boiling temperature. In the present cal- 
culations, the vaporization rate under normal evapo- 
rating conditions is evaluated by means of a simplified 
classical Z) 2 -theory: 


m k = 2nr k p g Df gs Sh ln(l + B k ) (64) 

where the Spalding mass transfer number, B k , is 
given by Eq. (24) and the Sherwood number, Sh, 
is given by 


Sh = 2(1 + 0.3 Re 1/2 Scl /3 ) (65) 

8.4 Internal Droplet Temperature Calculation 

Our experience with the validation studies 
showed us that there is a definite need to include 
a calculation involving the internal droplet tempera- 
ture valid under both superheat and normal evapo- 
rating conditions. In our present calculations, it was 
evaluated by means of a simple infinite conductivity 
model. 


dT k _ 3 [h,eff - h\ ds k 
dt 2 C p ir\ dt 

if T k < T b , and 


( 66 ) 


dT k 

dt 


3a., 


~{T k -T b ) 


r k piC p i 

( for the superheat model of [62, 63 — 64]) 


dT k _ 3A; dT 

dt r k piC p i dr 


(67) 


( for the superheat model of [74]) 


if T k > T b 


There is a need for understanding droplet va- 
porization behavior at supercritical conditions be- 
cause of the increasingly high operating pressures en- 
countered in some gas-turbine combustors. The en- 
gine operating pressures can exceed the supercriti- 
cal pressure of liquid fuels such as jet-a. Because 
of its practical and fundamental importance, several 
experimental and theoretical investigations were un- 
dertaken to understand droplet gasification occurring 
at supercritical conditions [19-20,82-86]. Most of the 
numerical investigations that appeared in the litera- 
ture to study the vaporization behavior of an isolated 
single spherical droplet were based on the coupled 
twophase, unsteady Navier-Stokes equations in 1-D 
spherical coordinates [20,83-86]. However, it is note- 
worthy that the resulting gas-liquid interface analysis 
at the droplet surface becomes extremely complicated 
for a multi-component mixture. Also, the numerical 
solution is complicated by several factors as it needs 
to take into account the high-pressure corrections as- 
sociated with various transport and thermo-physical 
properties. The supercritical droplet vaporization 
differs from low-pressure droplet vaporization models 
in several important ways. For example, the results 
of Zuo et al [85-86] showed that under low to mod- 
erate temperatures, first there would be an increase 
in droplet lifetime before it starts to fall off with an 
increase in pressure. Their results also showed that 
there would be a monotonic decrease in droplet life- 
time with increase in pressure at high ambient tem- 
peratures. 

However, none of these analyses could be di- 
rectly incorporated into a comprehensive spray calcu- 
lation because of any such detailed treatment would 
be prohibitively expensive in terms of the required 
computer CPU resources. Also, they fail to address 
how to handle the supercritical droplet vaporization 
in a real gas-turbine environment in which a droplet 
traverses through a non-stagnant gas. Since our main 
interest lies in implementing a viable high-pressure 
droplet vaporization model into our spray code, we 
attempted to incorporate some important aspects of 
high pressure modeling derived from Refs. [20,83-86] 
into the existing framework of our spray formulation. 
In what follows we describe some important aspects 
of supercritical droplet vaporization: 

(1) The non-idealities in gas-phase become in- 
creasingly more pronounced as the pressure ap- 
proaches a supercritical state. There are several 
widely used equations of state that provide accurate 
representation at high pressures, namely the Peng- 


N AS A/CR— 20 15-218918 


18 


Robinson (PR), Redlich-kwong (RK), and Soave- 
Redlich-Kwong (SRK). In our calculations, we im- 
plemented the PR EOS because of its simplicity and 
also known to provide accurate representation [85- 
86]. Also, it is known to provide reasonable results 
when calculating the phase equilibrium at the droplet 
surface. 

(2) At low pressure conditions, the solubility 
of the ambient gases in liquid phase can be ne- 
glected. The corresponding gas-phase composition 
at the droplet interface can be calculated by means 
of Raoult’s law, Eq. (42). At high pressures, the 
phase equilibrium calculations at the droplet surface 
need to take into account fugacity of each individual 
species. More on how to evaluate the phase quilib- 
rium at the droplet surface can be found in Section 
9.1. 

(3) Both liquid and gas compositions encoun- 
tered in a gas turbine combustion device are truly 
multi-component in nature. However, all the nu- 
merical investigations reported so far on supercritical 
droplet vaporization are based on a binary mixture 
[19-20,83-86]. Even for a binary mixture, the phase 
equilibrium calculations at high pressures present a 
formidable challenge. As it becomes very difficult 
to solve anything other than a binary mixture, we 
restrict the phase equilibrium calculations to a bi- 
nary mixture involving a combination of a single- 
component surrogate fuel and nitrogen. Such a com- 
bination seemed to provide a reasonable representa- 
tion for a fuel droplet vaporizing at high pressures 
[85-86]. 

(4) The transport properties in both liquid and 
gas phases become increasingly pressure dependent at 
high pressures. More on how we evaluate the trans- 
port properties at high pressures can be found in Sec- 
tions [4] & [7]. 

(5) As the droplet surface approaches a trans- 
critical state, the latent heat of vaporization dimin- 
ishes to zero. In our current calculations, the latent 
heat of vaporization is calculated by a formulation as 
described in Section 9.2. 

(6) The solubility of the ambient gases in liquid 
phase becomes increasingly important at high pres- 
sures. However in our present calculation, we ignore 
to take into consideration the multi-component na- 
ture of droplet behavior for the following reasons: 
For gas turbine combustors of our interest, the gas 
pressure seldom exceeds more than twice the criti- 
cal pressure of jet fuels. Within that pressure range, 
the liquid phase solubility of gas could be ignored 


since it could be shown from phase equilibrium cal- 
culations that the liquid-phase mass fraction of ni- 
trogen remains less than three percent even when the 
droplet surface temperature reaches near critical tem- 
perature. 

(7) Zhu et al. [85-86] studied the influence of 
gas-phase unsteadiness on droplet vaporization, and 
also quantified to some degree the resulting differ- 
ences from the quasi-steady and transient models. 
Their results show that unsteadiness seemed to per- 
sist over a wide region during a brief early transient 
period after a droplet is suddenly introduced into 
an otherwise stagnant gas. After this initial tran- 
sition period, some unsteadiness remains persistent 
in a small region closer to the droplet surface. Ini- 
tially, the quasi-steady model seemed to produce a 
smaller regression rate when compared with the tran- 
sient model. But in the later stages of droplet life- 
time, it seemed to produce a much higher regression 
rate. The influence of unsteadiness seemed to in- 
crease with an increase in both ambient pressure and 
temperature. However in our present calculations, 
we expect for the quasi-steady model to provide a 
useful approximation for the following reasons: (1) 
We are interested mainly in gas pressures not too far 
above the critical pressure of the liquid fuel, and (2) 
Also, we are interested in a droplet moving a strati- 
fied gas where the unsteadiness associated with initial 
transient can be neglected. It is because our droplet 
models are based on what happens after the initial 
atomization phase. So the neglected influence of un- 
steadiness is mostly limited to a small region near the 
droplet surface. 

9.1 Equilibrium Relations Valid at High Pressure 
Conditions 

At high pressures, the equilibrium needs to sat- 
isfy 

T l = T v (68) 

P L = P v (69) 

fi = fY => x i<t>i = Vi <t>Y ( 70 ) 

where fi is the fugacity of the ith species, Xi is the 
mole fraction species in the liquid phase, yi is the 
mole fraction species in the vapor phase, and 4>i is the 
corresponding fugacity coefficient of the ith species 
which is given by 
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A 


dZ 

dT 


In cf>i 


-(Z-l)-ln(Z-B*) + 
bm. 


2\/2 B* 


( 


bi 


—5i)ln 


Z + B*( l + y/2) 

Z + B*( 1 - y/2) 
(71) 


where = YJ&7K, ’ “ d 

Si = — 12i x j a y 2 ( 1 — kij)- There are two more 

additional equations that need to be satisfied at the 
droplet interface, 


T— 1 

II 

w- 

(72) 

£*< = 1 

(73) 


l 


For a given P and T of a binary mixture, there 
are six equations and six unknowns: Z v , Zi , yi, 2 / 2 , 
X\, and X 2 - The Peng- Robinson EOS yields solution 
for Z v and Z[ and the other four unknowns are cal- 
culated from the solution of Eqs. (70) to (73). The 
solution for this highly non-linear set of equations is 
obtained from the use of a Newton-Raphson iterative 
method. 


9.2 Calculation of Latent Heat of Vaporization at 
High Pressure 

The enthalpy of the species can be calculated 

from 

K^K-RT'f-^}^ (74) 

The related derivatives in the above equation 
are given by 


din 

DT 



da 

dT 


dA* 

dT 


lndZ_ _ (gf + + gf) 

b dT Z - B* Z 2 - B * 2 + 2 ZB* 

1 bi r w dA* , A%, Z + B*( l + y/2) 

2 y/2B* b dT T’ Z + B*{l-V2) 


dSi Z + B*(l + V2) 
-lr 


2V2B* dT Z + B*( 1-V2) 
Rfwi ( 0.45724T C j \ 1/2 


( 0.45724T ci ^ 1 

V Tp a ) 

EE-f— > 


= A* 


1 da 
adT 


|, 


(75) 

(76) 


(77) 

(78) 


3B* 

dT 


B 

T 


(79) 


(80) 


r 3A 

1 dT 


(Z - B*) + ^-[Z 2 - 2Z(1 + 3 B*) - A* + 2 B* + 3 B* 2 ] 
3Z 2 - 2(1 — B*)Z + A* — 2 B* - 3 B* 2 


(81) 


dSi 

~dT 


Si[~ 


1 da) /2 1 da 


!/2 dT 


1 da a X J 2 v— > 

+ 2 — — > 

adT ' a 


Vj ( i - k a)- 


da 


1/2 


dT 


(82) 


Finally, once the enthalpies in both phases are 
known, the latent heat of vaporization is simply cal- 
culated as follows 


Li =K- h\ (83) 

9.3 ncc_hp_eos_vle_table.dat: A Table Containing 
the Equilibrium Properties of a Binary Mixture 

Instead of performing these calculations on the 
fly, it would be a lot more economical to create a ta- 
ble containing the equilibrium properties of a binary 
mixture. Especially in gas-turbine combustor calcu- 
lations we need to tabulate these properties for one 
given pressure. Based on the calculations described 
in Sections 9.2 to 9.3, we developed the coding neces- 
sary to create a table, ncc_hp_eos_vle_table.dat, that 
contains values for the following variables: latent heat 
and mole fractions of both species, density, compress- 
ibility factor (Z), and specific volume in both liquid 
and gas phases. The table contains incremental val- 
ues over a wide range of expected liquid droplet tem- 
peratures for any given pressure. The values for any 
required temperature is calculated by means of linear 
interpolation. 

10 DETAILS OF ATOMIZATION 
MODELING 

Atomization refers to a process of the liquid jet 
breakup into droplets. There are many processes 
associated with the liquid jet breakup. They can 
broadly be classified into two regimes: (1) the in- 
ner nozzle flow, and (2) liquid behavior after exiting 
from the nozzle. In the inner nozzle flow, several 
factors (such as injector type, geometry, and size) in- 
fluence the conditions at the injector exit (such as 
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the velocity, the initial sheet or jet thickness, and 
the angle of droplet dispersion). One way to spec- 
ify the initial spray conditions at the injector exit is 
to rely on the widely-used correlations. For a bet- 
ter description, a more accurate analysis is needed 
which takes into account the physics associated with 
inside bubble growth, cavity formation, and internal 
turbulence. Once a liquid exits outside of the nozzle, 
it becomes unstable under the influence of aerody- 
namic instabilities and finally breaks up into droplets. 
The widely known aerodynamic instabilities are of 
Rayleigh-Taylor and Kelvin-Helmholtz kind [41, 43- 
44]. The Rayleigh-Taylor instability is due to inertia 
of the denser fluid opposing the system acceleration 
in a direction perpendicular to the interface of the 
denser fluid and the Kelvin-Helmholtz instability is 
caused by the viscous forces due to the relative mo- 
tion of the fluids [43-44]. When the maximum am- 
plitude of the most unstable wave reaches a critical 
value, some liquid is stripped of the main liquid core 
in the form of primary droplets. These droplets may 
further breakup into smaller droplets due to a process 
known as the secondary droplet breakup mechanism. 

Based on a linear instability analysis of a 2D 
viscous incompressible fluid moving thorough an in- 
viscid incompressible gas, Reitz and Bracco [41] char- 
acterized the break-up regimes to be four-fold: (1) 
Rayleigh breakup, (2) first wind-induced breakup, 
(3) second wind-induced breakup, (4) atomization. 
In the first two regimes, drops of sizes greater than 
or equal to the nozzle diameter are produced at dis- 
tances far from the nozzle exit. In the applications of 
our interest, the last two regimes are more important 
where droplets much smaller than the nozzle diam- 
eter are produced near the nozzle exit. The knowl- 
edge gained from the instability analysis of various 
kinds [41, 45-46] is combined with some experimen- 
tal observations to form the basis for various models 
developed for both primary atomization & secondary 
droplet breakup. In this aproach, the jet breakup 
is modeled first by making use of a drop represen- 
tation approach in which discrete parcels of liquid 
are injected in the form of blobs with a characteris- 
tic size representative of the nozzle diameter instead 
of tracking an actual intact liquid core that forms at 
the nozzle exit. In the case of a planar or conical 
liquid sheet, the discrete parcels essentially represent 
liquid ligaments. Before the jet breakup the discrete 
parcels stay inside of the liquid core or sheet but af- 
ter the jet breakup they move independently. The 
breakup criterion, atomization rate, drop size and ve- 


locity and the location of the newly formed droplets 
are primarily determined based on an instability anal- 
ysis derived from the conservation equations of mass, 
momentum and energy. The analysis of the jet or 
sheet breakup into ligaments or droplets, the strip- 
ping of the liquid into fragments or droplets, and the 
formation of more droplets from further breakup of 
ligaments or fragments are all described under the 
classification of primary atomization. 

Some of the large droplets that are formed im- 
mediately after the primary liquid jet breakup may 
further breakup into smaller droplets under the influ- 
ence of aerodynamic instabilities. The large droplets 
first tend to flatten under the influence of aerody- 
namic pressure. Then large amplitude long wave- 
length waves caused by drop deceleration induce 
a Rayleigh-Taylor instability on the flattened drop 
causing it to breakup further into several relatively 
large-size product droplets. While at the same time 
short surface waves induce a Kelvin-Helmholtz insta- 
bility on the windward side of the parent drop re- 
sulting in the generation of much smaller product 
droplets. The breakup of the larger droplets into 
smaller droplets is described under the classification 
of secondary droplet breakup. 

Most of the coding required for both the pri- 
mary atomization and secondary droplet breakup 
models was provided by CFDRC. These models are 
applicable under both non-superheat and superheat 
conditions. In the superheated regime, the thermo- 
physical liquid properties are replaced by the two- 
phase properties of a superheated fluid [74]. The CF- 
DRC atomization module contains the following four 
primary atomization models: (1) the sheet breakup 
primary atomization model, (2) the blob jet primary 
atomization model, (3) the air-blast atomization 
model & (4) the BLS (Boundary-Layer Stripping) 
primary atomization model, and it also contains the 
following three secondary droplet breakup models: 
(1) the Rayleigh-Taylor secondary droplet breakup 
model, (2) the TAB (Taylor Analogy Breakup) sec- 
ondary droplet breakup model, and (3) the ETAB 
(Enhanced TAB) secondary droplet breakup model. 
The choice between various models depends on the 
specific application. Complete details of the models 
contained in the atomization module can be found in 
the appendix. 

10.1 Flash Atomization Modeling 

There are three components to the atomization 
modeling of a superheat spray: (1) primary atomiza- 
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tion, (2) secondary droplet breakup, and (3) flash in- 
duced atomization. Further details of both primary 
atomization and secondary droplet breakup models 
are described in a previous section. For the super- 
heated sprays, we have to deal with a third compo- 
nent known as flash induced atomization. We have 
followed the modeling approach of Lee et al. [74] for 
the flash induced atomization. In what follows we 
describe some details of their modeling approach. 

Flash atomization is a result of rapid evapora- 
tion associated with strong nucleation and bubble 
growth occurring inside of a non-equilibrium two- 
phase fluid [74]. When the void fraction within the 
rapidly evaporating fluid reaches a certain critical 
void fraction, the inner liquid core would shatter into 
pieces. It seemed to have a value of 0.61 based on 
the research identified in [74], Incidentally, this con- 
dition for criticality also coincides with the packing 
limit of spheres beyond which it can not accommo- 
date any more bubbles and results in shattering [75]. 
In order to calculate this state of critical limit, they 
employed the homogeneous relaxation model (HRM) 
for calculating the mass fraction of gas, x , also known 
as quality. 


predicted by means of a surrogate model of Violi. 
Further details can be found in [74]. 

The flash induced atomization is combined with 
the primary and secondary droplet breakup models 
following the approach of Lee et al. [74]. The aero- 
dynamic instability mode will atomize the liquid core 
when the surface disturbance reaches a certain critical 
value and the flash induced atomization is responsi- 
ble for shattering the liquid core when the void frac- 
tion reaches the critical value. Therefore, the control 
mechanism depends on whichever mechanism is first 
responsible for reaching the critical condition. 

10.2 Some Effects of Superheat on Primary 
Atomization 

Flash evaporation is known to produce a fine 
spray with enhanced atomization, increase effective 
spray cone angle, and decrease spray penetration [61]. 
Here, we consider the effects of flash evaporation on 
the sheet breakup primary atomization model by fol- 
lowing the approach of Zuo et al [62]. Its effect on 
the initial droplet size generated immediately after 
the first ligament breakup, is given as a function 
of both engine pressure and a superheat parameter. 


f t = e = 0 o eV; (84) 

where the values for a (= —0.54), /3 (= —1.76), and 
0 o (= 6.51.10 04 s) are taken from [76], and the void 
fraction, e, and the degree of superheat, </>, are defined 
as follows: 


e = 


Pi- P . 
Pi Pv 


Pb{h ) - P 
Pc - Pb(h) 


(85) 


Once the critical limit for the void fraction is 
reached, the resulting droplets produced as a result 
of flash atomization are obtained by a droplet number 
density correlation developed by Kawano et al. [77]. 


n = Aem-Tbu b ) e at+b ( 86 ) 

From Eq. (86), the mean droplet size, R, is 
calculated by equating total liquid mass before and 
after the droplets are formed. 


R* = 


1 - e c 
(4/3)7 t 


(87) 


The solution for Eqs. (84) to (87) and all the 
JP8 fluid properties needed in this calculation were 
obtained by the computer coding provided by UTRC. 
Also, it is noteworthy that the properties for JP8 were 


= d i n(^-) 0 ' 27 [l-x(^) 0135 l 0<X(^)°' 135 <1 

(88) 

where di n is the corresponding droplet droplet size 
occurring under normal evaporating conditions with- 
out flash evaporation, and x, a superheat parameter, 
is defined as follows: 


I(T k ) - I(T b ) 

m) 


(89) 


where / is the internal energy of the liquid. Its value 
varies between 0 < x < 1 with x = 0 referring to zero 
flash evaporation and x = 1 to full flash evaporation. 

In Eq. (88), the increase in due to an in- 
crease in engine pressure by a factor of (p^— ) 0 ' 27 
is based on an experimental correlation obtained 
from Lefebvre [66]. It reflects the influence of cham- 
ber pressure on wave propagation as it damps wave 
growth. But the decrease by a factor of (1 — 
X(fk |m)°-i35) j s ( ] ue to a significant reduction in 
droplet size caused by both cavitation and bubble 
growth under flash evaporation conditions. It was 
introduced based on the experimental data obtained 
from VanDerWege et al [67] and Reitz [68]. 

As the liquid approaches boiling, it also causes a 
substantial decrease in both intact liquid core length 
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and core droplet size leading to a modification of the 
nominal cone angle, 0 , as given by 

e = 0 n + (144 — 0 n )x 2 (90) 

where 9 n is in degrees for a spray vaporizing under 
normal conditions without flash evaporation. This 
correlation was developed based on the experimental 
data of Reitz [68] . 

These two correlations were developed in con- 
junction with the sheet breakup primary atomization 
model but their application with other primary at- 
omization models need further investigation. 

11 DESCRIPTION OF INITIAL SPRAY 
CONDITIONS 

The spray computations facilitate the use of 
multiple fuel injectors. The same or a different type 
of liquid fuel can be specified for each one of differ- 
ent injectors. The initial droplet temperature is as- 
sumed to be the same for all different droplet groups 
of a given injector. The liquid fuel injection is simu- 
lated by introducing a number of discretized parcels 
of liquid mass at the beginning of every fuel- inject ion 
time step, A tu. The following three variables play 
an important role in simulating the injector initial 
conditions: 

• no_of_holes() - the number of holes per injector, 

• no_of_streams() - the number of droplet streams 
per hole - it is introduced to distinguish the 
initial velocity variation within different droplet 
classes arising from the geometric considerations 
of a chosen spray, & 

• no_o f .droplet .groups () - the number of droplet 
groups per stream - it is introduced to distin- 
guish the droplet-size variation within different 
droplet classes of a polydisperse spray. 

The total number of droplet groups introduced at 
the beginning of every different injection time step 
for a given injector is therefore equal to a value 
given by the multiplication of these three vari- 
ables. Further details on some of the spray in- 
put variables can be found in the spray input file, 
nccdnjector.in.l of Table 6, where the integer num- 
ber followed by the last dot represents the in- 
jector number, which in this case happens to be 
one. The table provides a complete description of 
the following input variables: nolc(n.i), out_string, 


(ymki(n_i,n_l), n_l=l,nolc(n_i)), tdrop(n_i), atomiza- 
tion(n.i), drop .breakup _model(n_i), spray_table(n_i), 
steady .spray .model, no .of .holes (nd) , 

no_of_streams(n_i), 

no_of_droplet_groups(n_i), lmdis(n_i), smdm(ni), 
cone(nd), size_min(n_i), size_max(n_i), stochas- 
tic(n_i), ((x_inj(n_i,nx), y_inj(n_i,nx), z_inj(n_i,nx), 
z_inj(n_i,nx), flowf(n_i,nx), v_inj(n_i,nx), 

alpha_inj(n_i,nx), beta_inj(n_i,nx), theta_inj(n_i,nx), 
dtheta_inj(n_i,nx), swlr_angle(n_i,nx)), 

nx=l,no_of_holes(n_i)), (atom_type(n_i,nx), 

breakup_type(n_i,nx), dia_hole(n_i,nx), 

delp_inj(n_i,nx), liq_vel(n _i,nx), gas_u(n_i,nx), 
gas_v(n _i,nx), gas_w(n_i,nx), pcl_start(n_i,nx), 
pcl_end(n_i,nx)), nx=l,no_of_holes(n_i)). 

The initial droplet distribution for a given in- 
jector could be specified by making use of one of the 
three available options: (1) by providing a complete 
specification of the initial conditions through a spray 
table, (2) by invoking certain pre-defined correlations, 
or (3) by chosing one of four available primary at- 
omization models. When Option (3) is selected, the 
logical variable, atomization(), is set to .true, or oth- 
erwise it is set to .false.. 

Option (1): 

When a spray table is defined, one should pro- 
vide the information on the (x,y z) components of 
the initial droplet location, the (u, v, w) compo- 
nents of the droplet velocity, and the initial mass flow 
rate associated with each different droplet group as 
described in ncc_spray .tabled of Table 7. This ta- 
ble provides a complete description of the following 
variables: nos(n_i), (ni,xx_inj, yy_inj, zz_inj, uu_inj, 
vv_inj, ww_inj, r.inj, fld.d), ,ni=l,nos(n_i)). The ini- 
tial inputs specified through a spray table should be 
representative of the integrated averages of the ex- 
perimental conditions [10-13]. 

Option (2): 

In this option, we need to specify several pa- 
rameters including no_of_holes(), no_of_streams(), & 
no_of_droplet_groups(). And depending on what is 
specified for the input of the integer variable, lmdis, in 
nccdnjector.in.l, the droplet-size distribution within 
each one of the streams is determined based on the 
following three choices: 

• In one choice, it is calculated based on a corre- 
lation typical of those widely used in describing 
the initial droplet size distribution [4]: 
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Table 6. Input file: ncc_liquid_injector.in.01. 

Input file content 

comments 

heading 

title of a description of this file 

heading 

title of property 

nolc(n_i) 

denotes the total number of initial liquid components 
in the nJ-th injector. 

heading 

title of property 

out -String 

chemical symbols of initial liquid species: e.g., C6H12. 

(ymki(n_i,n_l), 

n_l=l,nolc(n_i)) 

initial mass fractions of liquid species 

heading 

title of controlling parameters 

tdrop(n_i), 

isuperhd(nj), 

ijetapr(n_i), 

vmfdrop(n_i) 

tdrop(n_i) denotes the initial droplet temperature. 

If isuperhd(n_i) = 1, it invokes the superheat vaporization 
model of Zuo et al. [62], and Schmehl et al. [63-64]; If isuperhd(nj) 

= 2, it invokes the superheat vaporization model of Lee et al. [74]; & 
if isuperhd(n_i) = 0, it selects a non-superheat vaporization model. 

if ijetapr(nJ) = .true., all the JP8 fluid properties are calculated by 
the computer coding provided by UTRC as described in Lee et al. [74]. 

vmfdrop(n_i) denotes the initial quality ( = gas mass fraction) of a 
superheated fluid. It only applies when ijetapr(n-i) = .true. 

heading 

title of controlling parameters 

atomization(nd) , 
drop_breakup- 
model(n_i), 
line_injection(n_i) , 
spray _table(n_i), 
steady _spray_model 

If atomization(nJ) = .true., the initial droplet conditions are 
determined based on the use of a primary atomization model. 

Otherwise they are determined from either known 
experimental conditions or a widely-used correlation. 

If drop -breakup -model(n_i) = .true., it invokes the secondary 
droplet breakup option. 

if line_injection(n_i) = .true., it invokes a multi-point spray injection option. 

If spray -table = .true., initial droplet location, velocity, size, 
and mass flow rate are input through the file - ncc-liquid_table.in.l. 

Otherwise they are determined by some other means. 

If steady _spray .model = .true., it invokes a steady state spray model 
commonly used in many spray codes, where whenever the spray solver is 
called, after first introducing a new group of spray particles, it 
continues with the liquid-phase computation until after all the 
particles are taken out of the computational domain (note: It is NOT 
recommended to use this option as a steady state calculation could be better 
arrived at by making use of some features of the unsteady option). The solution 
from the unsteady option (steady-Spray_model = .false.) is determined based 
on the values assigned to the controlling time steps - dtml, dtil, 

& dtgl, which are internal to the spray solver. 
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Table 6. Input file: ncc Jiquid Jnjector.in.01 (continued). 

heading 

title of controlling parameters 

no_of_holes(n_i), 
no_of_streams(n_i), & 
no_of_droplet_groups(n_i) 

no_of_holes(n_i) = The number of holes per injector (note: 
when atomizationQ = .true, or spray .table = .true., the 
no_ofJioles(n_i) is set equal to 1). 


no_of_streams(n_i) = The number of droplet streams per hole. This 
variable is introduced mainly to distinguish the variation in the 
droplet groups based on angular orientation (note: when atomizationQ 
.true, or spray -table = .true., the no_of_streams is set equal to 1). 

no_of_droplct_groups(n_i) = The number of droplet groups per 
stream. This variable is introduced mainly to distinguish the 
droplet groups based on the droplet-size variation (note: 
when spray -table = .true, or atomizationQ = .true., all different 
injected droplets are lumped together into a single variable, 
no_oLdroplet -groups () ) . 

heading 

title of controlling parameters 

lmdis(n_i), 
smdm(nj), 
cone(nJ), 
size_min(n_i), 
size_max(n_i), & 
stochastic (n_i) 

If lmdis(n_i) = 1, it invokes the droplet size distribution as given 
by Eq. (91); If lmdis(n_i) = 2, it assumes a uniform droplet distribution 
between the maximum and minimum limits as specified by size_max() 

& size_min(); & If lmdis(n_i) = 3, it invokes the droplet size 
distribution as given by Eq. (92). 

smdm(n-i) = Sauter mean diameter 

If cone(nJ) = .true., it activates a 3D solid or hollow cone spray 
configuration as shown in Fig. 2. Otherwise it activates a 2D configuration 
depending on the value chosen for the logical variable - axisymmetric. If 
it is set equal to .true., it invokes the axis-of-symmetry case as shown 
in Fig. 3 otherwise it invokes the planar case as in Fig. 4. 

size_min(n_i) & size_max(n_i) = The variables are associated with 
the hndis(n_i) = 2 option. 

stochastic(n_i) is used in conjunction with Option (2) of Sec. 9. If set 
equal to .true., it introduces some randomness into the determination 
of initial droplet velocity distribution. 

heading 

title of controlling parameters 

((x_inj(n_i,nx), 

y-inj(n_i,nx), 

z_inj(n_i,nx), 

flowf(n_i,nx), 

vdnj (n T,nx) , 

alpha Jnj (n_i,nx) , 

beta_inj(n_i,nx), 

(x_inj(n_i,nx),y_inj(n_i,nx),z_inj(n_i,nx)) = spatial coordinates 
of the initial injector location. 

flowf(n_i,nx) = injector mass flow rate, (units - kgrn/s for 3D & 
axis-of-symmetry & kgm/s/m for 2D planar. REMEMBER FOR AXIS-OF- 
SYMMETRY, IT IS THE TOTAL FLOW RATE OVER 360 DEGS.) 
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Table 6. Input file: ncc_liquid_injector.in.01 (continued). 

thetadnj (n_i,nx) , 
dtheta Jnj (n_i,nx) , 
swlr_angle(n_i,nx)), 
nx=l,no_of_holes(n_i)) 

v_inj(n_i,nx) = droplet injection velocity, m/s. 
alpha_inj(n_i,nx) = angle of rotation from the x-y plane, 
beta Jnj (n_i,nx) = angle of rotation from the x-axis. 
theta_inj(n_i,nx) = cone angle. 

dtheta_inj(n_i,nx) = half-cone angle (note: Although dtheta_inj (n_i,nx) 
= thetadnj (n_i,nx)/2 for SOLID CONE SPRAY, it is invoked by 
setting dtheta dnj(n J,nx) either equal to 0 or thetadnj (nd,nx)/2). Refer 
to Figs. 2-4, for a better understanding of the angular representation. 

swlr_angle(nd,nx) = swirl angle. 

heading 

title of controlling parameters 

(atom_type(n_i,nx) , 

breakup _type(n_i,nx) , 

dia_hole(nd,nx) , 

delp_inj(n_i,nx), 

liq_vel(n_i,nx), 

gas_u(n_i,nx), 

gas_v(n_i,nx), 

gas w(n i,nx), 

pchstart(n_i,nx), 

pcl_end(n_i,nx)), 

nx=l,no_of_holes(n_i)) 

atom_type(nd,nx) = 0 => no primary atomization specified, 

= 1 => blob jet primary atomization model, 

= 2 => sheet breakup primary atomization model, 

= 3 => air blast primary atomization model, 

= 4 => BLS primary atomization model. 

breakup_type(nd,nx) = 0 => no secondary droplet breakup model specified, 
= 1 => Rayleigh-Taylor secondary droplet breakup model, 

= 2 => TAB secondary droplet breakup model, 

= 3 => ETAB secondary droplet breakup model. 

dia_hole(nd,nx) represents the injector orifice diameter (m). 

delp dnj(nd,nx) represents the pressure drop across the injector (Pa) 

(note: needed for pressure swirl atomization only). 

liq_vel(nd,nx) represents the liquid velocity from annulus (m/s) 

(note: needed for air-blast atomization only). 

gas_u(n_i,nx), gas_v(n_i,nx), gas_w(n_i,nx): gas velocity components 
in annulus adjacent to liquid film (m/s) (note: needed for air-blast 
atomization only). 

pcl_start(nJ,nx): starting parcel number for the annular air-blast injector. 
pcl_end(n_i,nx)): ending parcel number for the annular air-blast injector. 
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Table 6. Input file: ncc_liquid_injector.in.01 (continued). 

heading 

title of controlling parameters 

(line dnjection_type(n_i,nx), 

line_injection_type(n_i,nx) = 1 => line injection with full cone angle, 

xi_inj(n_i,nx), 

= 2 => line injection with half cone angle, 

yi_inj(n_i,nx), 

= 3 => circular injection with full cone angle, 

zLinj (n_i,nx) , 
xf_inj(n_i,nx), 

= 4 => circular injection with half cone angle. 

yf_inj(n_i,nx), 

(xi_inj(n_i,nx),yi_inj(n_i,nx),zi_inj(n_i,nx)) = In the case of line 

zLinj (n_i,nx) , 

injection, they represent the coordinates of the first end point of a line. 

ri_inj(n_i,nx)) 

Otherwise in the case of circular injection, they represent the 
coordinates of its center. 

In the case of line injection, (xf_inj(n_i,nx),yf_mj(n_i,nx),zf_inj(n_i,nx)) 
represent the coordinates of the second end point of a line. 

Otherwise in the case of circular injection, xf_inj (n_i,nx) & 
yf_inj (n_i,nx) represent the angular location of its two end points 
in degrees. 

In the case of circular injection, ri_inj(n_i,nx) represents its radius. 


Table 7. Input file: ncc_spray_table.in.01. 

Input file content 

comments 

nos(n_i) 

denotes the total number of droplet groups in the n_i-th 
injector. 

(ni,xx_inj,yy_inj, 

ni = number of the droplet group. Its value ranges between 

zz_inj ,uu_inj ,vv_inj , 
ww_inj ,r_inj ,fld_d) 

1 to nos(n_i). 

,ni=l,nos(n_i)) 

(xx_inj,yy_inj,zz_inj) = spatial coordinates. 
(uu_inj,vv_inj,ww_inj) = velocity components. 
r_inj = droplet size in radius. 

fld_d = mass flow rate of the ni-th droplet group (note: 
SUMMATION OF fld_d OVER ALL THE nos(n_i) 
DROPLET GROUPS IS EQUAL TO THE TOTAL MASS 
FLOW RATE OF THE INJECTOR). 
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Fig 2a. Geometrical details of fuel injection for 
a 3D solid or hollow cone spray. 



<|) = azimutal angle, deg 
number of rays = 8 (cj)=45) 

For a 3D cone, 
no_of_streams ( ) should 
defined as a multiple of 8. 
number along each ray = 
no_of_streams ( ) / 8 . 


Fig 2b. Initial spray particle orientation in a 
circular cross-section. 


N AS A/CR— 20 15-218918 


28 


0 = cone angle, deg 
d0 = half-cone angle 



Fig 3a. Geometrical details of fuel injection for 
an axis-of-symmetry case. 


• 00 

* For axis of symmetry: 

? no_of_streams ( ) is distributed 

• 01 between 01 to 00 . 


Fig 3b. Initial spray particle concentration in an 
axis-of-symmetry case. 


N AS A/CR— 20 15-218918 


29 


0 = cone angle, deg 



Fig 4a. Geometrical details of fuel injection for 
a 2D planar case. 


• 00 


♦ 01 For 2D: 


no 

of 

streams ( ' 

>/2 

is between 

01 

to 

00, and 



no 

of 

streams ( ' 

)/2 

is between 

LO 

to 

LI . 




* LI 

* 

* 


Fig 4b. Initial spray particle concentration in a 
2D planar case. 
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Figure 5 Droplet-size distribution. 
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where n is the total number of droplets and dn is 
the number of droplets in the size range between 
d and d + dd. This correlation also requires the 
specification of Sauter mean diameter, d^ 2 - Fig. 
5 shows the droplet size distribution generated 
by Eq. (91) for a case studied in [10]. The solid 
line shows the droplet number variation versus 
drop size and the dashed line shows the inte- 
grated mass variation with drop size. The drop 
size distribution within the spray is represented 
by a finite number of droplet classes as given by 
the variable, no_of_droplet_groups(). 


• In the second choice, the initial droplet sizes are 
distributed evenly between two specified maxi- 
mum and minimum drop sizes - size_max() & 
size_min() as specified in Table 6. And the size 
interval depends on the value specified for the 
variable, no_of_droplet_groups(). However, this 
option is applicable only in certain special cases. 


• In the third choice, the initial droplet size dis- 
tribution is calculated by making use of a cliped 
probability distribution function as given by 


P(X) = [ f(x)dx/ [ f{x)dx (92) 


to 


to 


where f(x) = (1.0 — exp(— x)(l+x+\x 2 + \x 3 ). 
After a random sampling of 0 < x < 12.0, the 
initial droplet diameter is determined as follows: 
(0.3 < 4 P(x) < 1. 5)^32- If the value for 4 P(x) 
goes beyond the limits of 0.3 and 1.5, the ran- 
dom sampling is reiterated until that number 
falls within the range of the lower and upper 
bounds. 

Depending on what is specified for the logical 
variable, cone, of Table 6, the droplet velocity dis- 
tribution amongst various streams of a given hole is 
calculated by assuming the spray to be either a solid 
or hollow cone spray. A graphical illustration of three 
different cone configurations are shown in Figs. 2 to 
4. Figs. 2a & 2b refer to a 3D case, Figs. 3a & 3b 
to an axisymmetric case, and Figs. 4a & 4b to a 2D 
planar case. It is noteworthy that in an axisymmet- 
ric case, the x-axis is assumed to be aligned with the 
axis-of-symmetry. 

Fig. 2a shows the geometric details of a hollow 
cone spray in 3D where 6 is known as the cone angle, 
dO is the half-cone angle, and for a solid cone spray 
dO = d/2. And a represents the angular rotation 
from the x-y plane and /3 is the angular rotation from 
the x-axis. Fig. 2b shows initial spray stream orien- 
tation in a circular cross section. We try to simulate 
the spray by a finite number of streams as given by 
the variable, no_of_streams(). Each one of the dark 
circles in Fig. 2b represents a different stream. The 
streams are distributed evenly along each one of the 
different rays which are separated from each other 
with an angle of separation as given by (j>. Currently, 
we have hard-coded the angle of separation to be 45° 
which means we have restricted the number of rays 
to be eight in 3D. Therefore, when specifying a num- 
ber for no_of_streams(), it should be borne in mind 
that it should be a multiple of eight. And the num- 
ber of streams along each one of the rays, therefore, 
becomes equal to no ofstreams()/8. In Fig. 2b, the 
no_of_streams() has a value of 32 and, therefore, the 
number of streams along each one of the rays for this 
case is 4. 

Fig. 3a. shows the geometric details for an 
axisymmetric case. Since the computations are per- 
formed only in the first quadrant of the x-y plane, all 
of the specified number of streams, no_of_streams(), 
are distributed evenly between the lower and upper 
halfs as shown in Fig. 3a and 3b. Here, the upper half 
refers to the region between OI to OO and the lower 
half refers to LO to LI. It is noteworthy that each 
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droplet group in a given stream represents a circular 
ring of liquid for an axisymmetric case. 

The geometric details for a 2D case are shown 
in Figs. 4a & 4b. Here, it is noteworthy that 
each droplet group in a given stream represents a 
planar sheet of liquid. All the specified number of 
streams, no_of_streams(), are distributed evenly over 
both sides of the cone center as shown in Fig. 4b. 

For each one of the different injector holes, 
no_of_holes(), it requires the specification of the fol- 
lowing parameters as described in ncc injector. in. 1 
of Table 6 (However, because of the geometric dif- 
ferences that exist between 3D, 2D, and axisymmet- 
ric configurations, some of the input parameters may 
have different units as noted below): 

• The initial (x, y, z) coordinates of the hole loca- 
tion. 

• The mass flow rate per hole - however, the defi- 
nition of the units for the injector mass flow rate 
per hole differs: it is kgm/s for 3D & the axis- 
of-symmetry and kgm/s/m for 2D planar (note: 
In an axis-of-symmetry, the specified mass flow 
rate per hole refers to the entire mass flow rate 
over 360 degrees). 

• The following variables define the angular ori- 
entation: otinj =angle of rotation from the x-y 
plane, j3i n j = angle of rotation from the x-axis, 
0 in j = cone angle, & d9i n j = half-cone angle 
(note: Although d,0i n j = 9i n j / 2 for a solid cone 
spray, a specified value of zero for d9i n j also in- 
vokes a solid-cone spray configuration). 

• The variable, swlr_angle(), allows a means to 
specify the tangential component in the case of 
both 3D and axisymmetric sprays. 

Stochastic Option 

In the above description of Option (2), the ini- 
tial droplet conditions are determined based on a de- 
terministic approach. However, such an approach 
when employed in 3D may lead to time consuming 
calculations requiring the use of many droplet groups. 
With this in mind, we wanted to introduce some ran- 
domness into the calculation of the initial droplet ve- 
locity distribution as an option. This option can be 
invoked by setting the logical variable, stochasticQ, 
to be .true.. In a stochastic approach, the initial an- 
gular orientation of a particle could be randomized in 
a number of different ways. The obvious way would 


be to determine the initial droplet direction based on 
a complete randomization in a field as defined by the 
geometric parameters: solid cone angle, 9, half cone 
angle, d9 and/or azimuthal angle, 0. Another way 
is try to retain the basic features of the determinis- 
tic approach where the initial droplet directions are 
taken as the means in a random field with the ran- 
domness localized to within the surrounding angular 
sub-segment of the respective droplet as represented 
by 59, 5d9, 5950, or 5d950. We have chosen to go 
with the later option as it preserves the same basic 
structure but provides some benefits of the stochastic 
approach. 

Option (3): 

Here many aspects of fuel injection remain the 
same as in Option (2) but the initial droplet size 
and velocity distribution may differ depending upon 
the primary atomization model selected. More de- 
tails on the available primary atomization models 
can be found in the appendix. For each one of the 
different injector holes, no_of_holes() it requires the 
specification of all the parameters as in Option (2). 
However, it doesn’t make use of the values set for 
size_max() & size_min(). Also, it requires the spec- 
ification of some additional parameters as described 
in ncc Jnjector.in.l of Table 6. However, there exist 
some differences with regard to their usage depending 
upon the primary atomization model. Some major 
differences are highlighted below: 

• The integer variable, atom_type() of Table 6, is 
used to select the primary atomization model of 
one’s choice. 

• dia_hole() needs to be specified for all primary 
atomization models. 

• delp_inj() needs to be specified only with the 
sheet breakup primary atomization model and 
it is used in Eq. (9) of the appendix. 

• liq_vel(n_i,nx), gas_u(n_i,nx), 

gas_v(nJ,nx), gas_w(nJ,nx), pcLstart(n_i,nx), & 
pcLend(n_i,nx)) are required only for the air- 
blast primary atomization model. 

• In all of the primary atomization models except 
for the air-blast, the total number of droplet 
groups (particles) to be injected at the time 
of every injection period is specified though 
the use of the variable, no_of_droplet_groups(). 
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Fig. 6 A graphic illustration of a multi-point (line) injection. 



Fig. 7 A graphic illustration of a multi-point (circular) injection. 
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For the air-blast, it is given by pcl_end(n_i,nx)- 
pcLstart (n_i,nx) + 1 . 

• In this option, all of the injected droplet groups 
as specified by either no_of_droplet_groups() 
or pcLend(n_i,nx)-pcLstart(n_i,nx)+l are dis- 
tributed randomly across the entire cross-section 
of any selected 2D, axisymmetric, or 3D spray 
configuration. So it is recommended that 
no_of_streams() be set equal to one as no dis- 
tinction is made between streams. 

The secondary droplet breakup option can be 
invoked regardless of how the initial droplet condi- 
tions are generated using Options (1) to (3). The sec- 
ondary droplet breakup option can be invoked by set- 
ting the logical variable, drop .breakup _model(), to be 
.true.. And then the integer variable, breakup _type() 
of Table 6, is used to select the secondary droplet 
breakup model of one’s choice. More details on the 
available secondary droplet breakup models can be 
found in the appendix. 

12 MULTI-POINT (LINE OR CIRCULAR) 
SPRAY INJECTION 

In order to perform some reacting spray calcu- 
lations involving a GE Taps concept which was devel- 
oped as a part of the NASA’s fundamental aeronau- 
tics/supersonic initiative on high altitude emissions, 
we implemented some new methods of spray injection 
into our spray code: (1) a line injection, and (2) an 
injection along a circular arc. 

In line injection, a new droplet is injected ran- 
domly between two end points of a line as specified 
by (Xi,inj and (& f ^inj ; V f ,inj 5 Z f ^inj') ■ The 

initial droplet velocity is defined by a velocity magni- 
tude as specified by Vi n j and its orientation as spec- 
ified by an angle, fii n j , in degrees. Within the line 
injection, we have added two different options: (1) 
line injection with a full cone angle, and (2) line in- 
jection with a half cone angle. In both these options 
the angle is given by 9i n j in degrees, and a newly in- 
jected droplet is chosen randomly to fall within the 
specified cone angle. A graphic illustration of the line 
injection with a full cone angle is shown in Fig. 6. 

In circular injection, a new droplet is injected 
randomly along a circular arc between two end 
points along the circumference of a circle. For this 
case, the coordinates of the circular arc are speci- 
fied by the following parameters: its center location, 
( Xi,inj,yi,inj’ z i,i n j ), its radius, r i>inj , and the angu- 
lar location of its two end points, Xf^ n j and y/,mj in 



Fig. 8 A vector illustration used in the particle search 
analysis. 

degrees. The initial droplet velocity is defined by a 
velocity magnitude as specified by Vi n j and its ori- 
entation as specified by an angle, /3i n j, in degrees. 
Within the circular injection, we have added two dif- 
ferent options: (1) circular injection with a full cone 
angle, and (2) circular injection with a half cone an- 
gle. In both these options the angle is given by 9i n j 
in degrees, and a newly injected droplet is chosen 
randomly to fall within the specified cone angle. A 
graphic illustration of the circular injection with a 
full cone angle is shown in Fig. 7. 

It is also important to note that the initial 
droplet size is determined based on a given value for 
SMD and a clipped probability density function as 
given by Eq. (92). And the number of droplets in 
a given parcel, (rifc), is determined based on the ini- 
tial liquid mass flow rate and the specified value for 
no_of_droplet_groups(n_i). 

13 SPRAY SOLUTION ALGORITHM 

In order to evaluate the initial conditions that 
are needed in the integration of the liquid-phase equa- 
tions, we first need to know the surrounding gas- 
phase properties at each particle location. But in or- 
der to evaluate the gas-phase properties, it is first nec- 
essary to identify the computational cell in which a 
given particle is located. It is a trivial task to track a 
particle in the regular rectangular coordinates. How- 
ever, the particle tracking becomes complicated when 
the computational cell is no longer rectangular in the 
physical domain and it becomes even more compli- 
cated when the particle search is performed within 
the context of parallel computing. 

We have developed and implemented an effi- 
cient particle-tracking algorithm for use with parallel 
computing in an unstructured grid. The search is ini- 
tiated in the form of a local search originating from 
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the computational cell in which the same particle was 
found to be located in the previous time step. The 
location of the computational cell is then determined 
by first evaluating the dot product of x pc . a n = \x pc \ 
|a„| cos ((/)), where x pc is the vector defined by the 
distance between the particle location and the center 
of the n-face of the computational cell, a n is the out- 
ward area normal of the n-face as shown in Fig. 8, 
and (j) is the angle between the two vectors. 

A simple test for the particle location requires 
that the dot product be negative over each and ev- 
ery one of the n-faces of the computational cell. If 
the test fails, the particle search is carried over to the 
adjacent cells of those faces for which the dot prod- 
uct turns out to be positive. Some of those n-faces 
might represent the boundaries of the computational 
domain while the others represent the interfaces be- 
tween two adjoining interior cells. The search is first 
carried over to the adjacent interior cells in the di- 
rection pointed out by the positive sign of the dot 
products. The boundary conditions are only imple- 
mented after making sure that all other remaining 
possibilities point towards a search exterior of the 
computational domain. This implementation ensures 
against any inadvertent application of the boundary 
conditions before the correct interior cell could be 
identified. 

After the gas-phase properties at the particle 
location are known, the solution for the ordinary dif- 
ferential equations of particle position, size, and ve- 
locity are advanced by making use of a second-order 
accurate Runge-Kutta method. The partial differen- 
tial equations governing the droplet internal thermal 
and mass transport are integrated by making use of 
a fully implicit Newton-Raphson iteration method. 

Finally, the liquid-phase source terms of the gas- 
phase conservation equations (1-4) are evaluated by 
making use of a time-averaging method. 

14 DESCRIPTION OF COMPUTER 
FLOWCHART, & THE TIME-AVERAGING 
SCHEME USED IN THE GAS-PHASE 

SOURCE-TERM CALCULATIONS 

• In order to know more about the time-averaging 
method, we need to know first about the three 
different time steps that are internal to the spray 
code: A t mU A t iU and A t gi . 

A t m i - the actual time step used in integrating 
the liquid-phase equations which is determined 
based on the smallest of the different time scales 


associated with various rate-controlling phenom- 
ena of a rapidly vaporizing droplet. 

A tn ~ the injection time step. It is the time 
step at which a new discretized parcel of different 
droplet groups are introduced into the computa- 
tion. 

A t g i - the global time step. Its introduction 
seems to provide better convergence in both un- 
steady and steady-state computations. 

• When the spray solver is called it advances the 
liquid phase equations over a number of itera- 
tions as determined by the ratio of At g i/ At m i. 

• It then evaluates the time-averaged contribution 
of the liquid-phase source terms, S g i, of the gas- 
phase governing equations (1-4) as follows: 

s„ = E w s " < 93 > 

LALqi 
m—1 u 

where 

M 

Y A t ml = A t gl (94) 

m—1 

• The values for A t m i, A t g i, & Atu are specified 
in the input file, ncc_liquid_solver.in, of Table 8. 

• In steady-state computations, it is recommended 
to use for both A t g i and Atu a value of about 
1 ms which is roughly equivalent to the aver- 
age lifetime of the droplets for a typical reacting 
spray encountered in conventional low-pressure 
gas-turbine combustors. 

The averaging scheme could be explained better 
through the use of a flow chart shown in Fig. 9. 
The main spray solver is invoked with a controlling 
routine, DCLR, which, then, executes the following 
steps: 

1. It first initializes the source terms to zero. 

2. Updates the global time, t g i, based on At g i. 

3. Checks to see if t m i < t g i < t m i + A t m i. If it 
is, it returns control over to the calling routine 
and supplies the other flow solvers, e.g., flow or 
EUPDF, with the source terms, S g i, of Eqs. (1)- 
(4). If not, it proceeds with the next step. 
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Table 8. Input file: nccdiquidjsolver.in. 

Input file content 

comments 

heading 

title of controlling parameters 

ldread,ispray_mod, 
(when_start_spray(n_i), 
n_i= 1 ,no_of -injectors) 

If ldread = .true., restarts the calculation from the data stored from 
the previous computation. Otherwise initiates a new spray computation. 

ispray_mod controls the calls to the spray solver. The spray solver is called 
once at every other CFD iteration as given by the number, ispray_mod. 

when_start_spray() represents the CFD iteration where the computations for 
the nJth injector are initiated. 

heading 

title of controlling parameters 

iswim, 

walLfilm_thickness, 
offset, & 

coef_of_restitution 

iswim = 0, default spray wall interaction model based on a single outcome of 
droplet moving along a wall, & = 1, new model based on four different 
possible outcomes. 

wall_film_thickness = 0, if dry wall, any positive number if wet wall, 
offset = offset fraction for boundaries. 

coef_of_restitution = coefficient of restitution at impact wall for drop. 

heading 

title of controlling parameters 

hp.eosl 

hp_eosl = 0, default (low pressure) spray model, & 

hp_eosl = 1, high pressure spray model as described in Sec. 9. This 
option requires the data contained in the file, ncc_hp_eos_vle_table.dat, 
as described in Section 9.3. 

heading 

title of controlling parameters 

dtml, dtgl, & 
dtil 

dtml = time step for advancing the liquid phase equations. 

dtgl = global time step. Whenever the spray solver is called, it advances 
the spray computations over a period of dtgl before returning control 
over to the calling routine. To be more precise, it advances the 
liquid phase equations over a number of time steps as determined by 
(dtgl /dtml). 

dtil = injection time step. It is where a new group of droplets are 
introduced into the computation. 
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Fig. 9 The flowchart of LSPRAY-V. 
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4. Checks to see if it it is time to introduce a new 
group of particles. 

5. Proceeds with solving the liquid-phase equations 
with calls to the following routines: 


It can be determined based on the the known initial 
particle location, x ik , particle velocity, V_ ik , and other 
geometrical considerations of the grid cell. Based on 
vector analysis, we have determined it as follows 


• Calls the particle tracking routine and as- 
signs particles based on the parallel strategy 
implemented. 

• Interpolates gas-phase properties at the 
particle location. 

• Advances liquid-phase equations and, then, 
deletes any particles that are no longer 
needed in the computations. 

6. Evaluates the liquid-phase source-term contribu- 
tions, S m i, of Eq. 93. 

7. Updates the time, t m i, based on A t m i. 

8. It then goes back to step (3) and repeats 
the whole process again until the computa- 
tions are completed over a global time step: 

Em=l = A t gl . 

15 IMPLEMENTATION OF BOUNDARY 
CONDITIONS 

The spray code supports all of the boundary 
conditions that exist in the current version of the 
NCC CFD module. First, we are going to provide 
a summary of the droplet and wall interaction mod- 
els. There are several possible outcomes following a 
droplet impact with a wall. They can be broadly clas- 
sified into four categories: (1) the droplet may move 
along a wall after its impact but still keeps on vaporiz- 
ing, (2) the droplet may rebound after its impact with 
the wall, (3) the droplet may stick to the wall after 
its impact but keeps on vaporizing, or (4) the droplet 
may be shattered after its impact with the wall. The 
first outcome is the default boundary condition that 
was originally implemented into our code. It is based 
on the assumption that a droplet, after having lost 
most of its momentum upon impingement with the 
walls, moves along the wall surface with a velocity 
closer to that of the surrounding gas-film. The last 
three outcomes are based on the models described in 
Refs. [80-81]. They are implemented partially based 
on the coding received from CFDRC. 

In order to implement the droplet and wall in- 
teraction models, there is a need to determine the in- 
tersection location of a particle crossing a wall, x bl k . 


—bl k — (*Ufc 4“ ^ Uinki Uik 4“ ^ V ink: %ik 4“ U Wi nk ) 

(95) 

where 


* Qin • i—cwf — ik ) 

u = V ’ 

—n ’—ink 

Y-ink > is the velocity normal (= V ik /\V ik \), and x cwf 
is the center location of the wall face of the grid 
cell. In what follows, we define various variables - 
Reynolds number, Re SWik , Weber number, We alU:k , 
Ohnesorge number, Oh aW}k: surface energy, E aw , and 
viscous dissipation, Vdi S - used in the spray-wall in- 


raction modeling. 



R&sWik — 

2pfcU.|UJ 

Pk 

(96) 

&sw,k — 

2pkfk\y_ik\ 2 

& k 

(97) 

Oh S w,k — 

\J~W &sw,k 
R^sw,k 

(98) 

E sw = no k (2r k /3max) 2 (,l - cos a)/4 

(99) 


where 


^JvTY-VYe^k 4We SWjk 

Pmax — \ ' 

3(1 - cos a) yj Re SWi k 

the droplet impact angle, a, is given by 


a = — — cos 


7, 


and the droplet frequency, 7, is 

— x ik ).a n 


7 = 


2 r k 


and 


Ydi a 


2nr k pk\Yi k \ 2 {2r k l3 

max ) 


3 y/R 


&sw,k 


(100) 


The droplet outcome after the interaction with 
the wall is determined based on the following criteria: 
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Shattering Droplet: 

A droplet shatters upon impact with a wall 
when We SWik > W e cr u, where the critical Weber 
number is defined as follows 


We crit = 9.9 10 Oh 2 8 

(101) 

The average size of the shatterred droplet 
impact is given by 

size after 

T ave,shat — 'WlQ'3'(j'll,shat 1 ^ 12 , shat) 

(102) 

where rn s h a t and ri 2 tS hat represent the two limits of 

droplet sizes as determined by 



(103) 

Z 1 ,s hat r-v | T r io 

2 Pk\y_ik\ 

4nr 3 k a k 

' 12, shat — j-, 

SU) 

(104) 

Once r ave s h a t is determined, the total number of 

shattered droplets produced is given by 


N s hat - int[( ) 3 ] 

^ ave,shat 

(105) 


However, it is important to note that the number of 
shattered droplets generated is limited to 2 as shown 
below 


N shat = min ( 2 , int [{ — — — ) 3 ]) (106) 

ave,shat 

The actual size of the shattered droplets is calculated 
by introduing some randomness as shown below 

T shat — T’ave,shat (0.85 + 0.30 RND) (107) 

where RND is a random number between (0,1). The 
droplet position and velocity of the shattered droplets 
is determined by the following criteria. 

Upon impact with a wall, a shattered droplet 
emerges from the wall in a direction determined by 
the modified reflection normal vector, x rnk . 

—rnk 2(0.5 RND)(x ink Cresix.ink'—n)—n) 

(108) 

where C res stands for coefficient of restitution. Its 
value depends on the properties of the wall but has a 
value of 2 under normal reflecting conditions. Also, 


RND is added to introduce some randomness into 
the direction of reflection vector. Once x rnk is de- 
termined, the new droplet position and velocity are 
determined. 

If Xrnk'—n ^ 0 ; 


i— nk 

\Y-ik l^rn/c 

(109) 

+ 1 %bl,k ~~ 

^.ik\^-rnk ^of fQin 

(110) 


where C Q f f stands for offset and is usually assigned a 
tiny value but in our present calculations it is assigned 
a value of 0.01 \x cw f — x cgc \, where x cgc is the center 
location of grid cell. 

but if x rnk .a n < 0, 

v nk = io - 06 « n ( 111 ) 

x nk = x b i k +C 0 ffa n (H2) 

Once the droplet size, location, and velocity of 
the shattered droplets are determined, the rest of the 
droplet properties can easily be deduced from the ini- 
tial conditions of incident droplet. 

Sticking Droplet to the Wall: 

A droplet sticks to the wall when E sw > Vdis- 
After being stuck to the wall, the impacted droplet 
size remains the same but continue to vaporize in 
time. The droplet velocity and position are assigned 
the following values: 


v nk = io - 06 « n 

(113) 

—nk = %bl,k + Coffin 

(114) 


The rest of the droplet properties are deduced 
from the initial conditions of incident droplet. 

Rebounding Droplet: 

If none of the above conditions are met for the 
outcomes of either droplet shattering or droplet stick- 
ing to a wall, the droplet is assumed to rebound after 
its impact with a wall. After the rebound, the droplet 
size remains intact as in the incident droplet. Both 
the droplet velocity and position are calculated as 
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in the droplet shattering but randomness is removed 
from the reflection vector. So it is calculated as in 
Eqs. (109-112) but RND is assigned a value of zero. 
The rest of the droplet properties are deduced from 
the initial conditions of incident droplet. 

Moving Droplet Along the Wall: 

For those droplets that were produced as a re- 
sult of droplet shattering, they are assumed to move 
along next to the wall surface after its subsequent 
impact with a wall. Such a droplet would experience 
no change in its size after its impact but it is allowed 
to keep on vaporizing as it moves along the wall sur- 
face. The droplet velocity and position are calculated 
as follows 

V nk = V gs (115) 

—nk = %U,k + °offa n (116) 

where V_ gs is the velocity of the the surrounding gas. 

Finally, in what follows we highlight how some 
other boundary conditions are implemented: 

• The implementation of the periodic boundary 
conditions becomes rather complicated. It is be- 
cause one needs to keep track of a particle leav- 
ing the computational domain from one periodic 
boundary, and for every particle that leaves the 
domain, a second particle reenters the compu- 
tational domain through a corresponding second 
periodic boundary. Also, one needs to take into 
account the possibility of the two computational 
cells where the particles leave and reenter the 
domain being assigned to a different processor. 
We incorporated the periodic implementation as 
a part of the particle search algorithm. It is im- 
portant to note that the boundary conditions are 
implemented with the help of some appropriately 
defined transformation matrices. Also, the par- 
ticle velocities of the entering particle need to be 
adjusted accordingly based on the orientation of 
respective periodic boundaries. 

• The symmetric boundary condition is imple- 
mented in such a way to satisfy the criterion that 
for every particle crossing the symmetry line, a 
similar one re-enters the domain in a direction 
determined by the reflection vector. 


• When the particles move out of the exit bound- 
ary, they are taken out of the computation. 


16 DETAILS OF COUPLING BETWEEN 
LSPRAY-V AND OTHER SOLVERS 

• The spray code is designed to be a stand-alone 
module such that it could easily be coupled with 
any other unstructured-grid CFD solver and the 
same holds true for EUPDF. (However, some 
grid-related parameters such as area vectors, grid 
connectivity parameters, etc. need to be sup- 
plied separately). 

• The spray solver needs information on the gas 
velocity and scalar fields from the other solvers 
and, then, it in turn supplies the liquid-phase 
source terms. 

• The PDF solver needs information on the mean 
gas velocity, turbulent diffusivity and frequency 
from the CFD solver and the liquid-phase source 
terms from the spray solver, and then it in turn 
provides the solution for the scalar (species and 
energy) fields to the flow and spray solvers. 

• It should also be noted that both the PDF and 
spray solvers are called once at every other spec- 
ified number of CFD iterations. 

• All of the three solvers (LSPRAY-V, EUPDF, 
and CFD) are advanced sequentially in an iter- 
ative manner until a converged solution is ob- 
tained. 

• All three codes (EUPDF, CFD, and LSPRAY- 
V) were coupled and parallelized in such a way 
to achieve maximum efficiency. 

The coupling issues could be better understood 
through the use of a flow chart shown in Fig. 10. 
It shows the overall flow structure of the combined 
CFD, LSPRAY-V, and EUPDF modules. Both the 
PDF and spray codes are loosely coupled with the 
CFD code. The spray code is designed in such a 
way that only a minimal amount of effort is needed 
for its coupling with the flow and PDF solvers. The 
present version of the spray module relies entirely on 
the use of Fortran common blocks for its informa- 
tion exchange with other modules. Even this reliance 
should entail only few changes to be made within the 
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Fig. 10 The overall flowchart of the combined CFD, LSPRAY, and EUPDF solvers. 
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spray code for its linkage with different solvers. The 
PDF code is also structured along similar coupling 
principles. 

The flow chart of Fig. 10 contains several blocks 
- some shown in black and/or solid lines and the oth- 
ers in color and/or dashed lines. The ones in solid 
blocks represent the flow chart that is typical of a 
CFD solver. The ones in dashed blocks represent the 
additions arising from the coupling of the spray and 
PDF solvers. 

The coupling starts with the calling of the sub- 
routine - spray _pdf_read_parameters, which then 
reads the spray control parameters from the input hie, 
nccJiquidsolver .in of Table 8. This table provides 
a detailed description of the following input hie vari- 
ables - ldread, ispray_mod, (when _st art _spr ay (n _i ) , 
n_i=l,no_of_injectors), iswim, wall_hlm_thickness, 
offset, coef_of_restitution, hp eosl, dtml, dtgl, & 
dtil. The coupling is then followed by the call- 
ing of the pdf_int_rerun subroutine. It initial- 
izes PDF computations and also it may restart 
the PDF computations if needed from the data 
stored from a previous iteration. Similarly, we 
call spray int rerun for the spray computations. 
It is noteworthy that the spray computations can 
be restarted by reading the data from the restart 
hies: nccJ.iquidjparams.out , nccJiquidjresults.db , 

& nccdiquidjresultsjmc.db. The restart capability 
is invoked by setting the logical variable, ldread, of 
the input hie, nccJiquidsolver.in of Table 8, to be 
true. Otherwise, the spray computations are initial- 
ized to start from the beginning. Then, the coupling 
proceeds with the calling of the following subroutines: 
dclr for integrating the spray calculations and eupdf 
for the Monte Carlo PDF. The input variable, is- 
pray_mod of Table 8, controls the calls to the spray 
integrating routine. The spray solver is called once at 
every other number of CFD iterations as specihed by 
ispray_mocl. And the hrst call to the spray solver is 
controlled by the input variable, when start spray () 
of Table 8, which represents the starting CFD itera- 
tion number from which the spray computations are 
initiated. Finally, the coupling ends with the call- 
ing of a subroutine, spray _pdf_output, which will 
create a set of new restart files. 

17 DETAILS OF FORTRAN 
SUBROUTINES & FUNCTIONS 

Table 9 provides a list of all the Fortran sub- 
routines developed as a part of the spray module. 


This table also provides information on all the For- 
tran functions. It also describes the purpose of all 
the individual subroutines and functions. Similary, 
Table 10 provides a list of all the Fortran functions 
and subroutines associated with the implementation 
of the high pressure equation of state as described in 
Sec. 3, and the high pressure behavior on the gas- 
phase transport properties as in Sec. 4. 

18 DETAILS OF PARALLELIZATION 

There are several issues associated with the par- 
allelization of both the spray & PDF computations. 
The goal of the parallel implementation is to extract 
maximum parallelism so as to minimize the execu- 
tion time for a given application on a specified num- 
ber of processors [37] . Several types of overhead costs 
are associated with parallel implementation which in- 
clude data dependency, communication, load imbal- 
ance, arithmetic, and memory overheads. The term 
arithmetic overhead is the extra arithmetic opera- 
tions required by the parallel implementation. Mem- 
ory overhead refers to the extra memory needed. Ex- 
cessive memory overhead reduces the size of a prob- 
lem that can be run on a given system and the 
other overheads result in performance degradation 
[37]. Any given application usually consists of sev- 
eral different phases that must be performed in cer- 
tain sequential order. The degree of parallelism and 
data dependencies associated with each of the sub- 
tasks can vary widely [37]. The goal is to achieve 
maximum efficiency with a reasonable programming 
effort [37]. 

In our earlier work, we discussed the parallel 
implementation of a spray algorithm developed for 
the structured grid calculations on a Cray T3D [10]. 
These computations were performed in conjunction 
the Monte Carlo PDF method. The parallel algo- 
rithm made use of the shared memory constructs ex- 
clusive to Cray MPP (Massively Parallel Processing) 
Fortran and the computations showed a reasonable 
degree of parallel performance when they were per- 
formed on a NASA LeRC Cray T3D with the number 
of processors ranging between 8 to 32 [10]. Later on, 
the extension of this method to unstructured grids 
and parallel computing written in Fortran 77 with 
PVM or MPI calls was reported in [11-15]. The lat- 
est version in Fortran 77/90 offers greater computer 
platform independence. In this section, we only high- 
light some important aspects of parallelization from 
Refs. [10-15], 
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Table 9. Description of LSPRAY-V Fortran subroutines & functions. 


Function 

Purpose of the Function 

blasiu(x) 

This function returns a solution for the function, /(!?&), 
of Eq. (21) for use in computing the droplet regression rate. 

Subroutine 

Purpose of the Subroutines 

axisyrmpart _adj ust 

This routine was written by Jeff Moder and essentially 
performs the same functions as the intpla routine but 
modified for application with a 2d-axisymmetric case. 

calc_F_suvw() 

It evaluates the right-hand sides of ODEs for both 
droplet size and velocities, and then provides the 
solution to the calling routine, chaslv. 

chaslv 

This routine has two main functions: 

(1) It integrates the liquid phase equations. 

(2) It removes the particles that are no longer needed 
in the computation. 

dclr 

This routine is called once at every other CFD iteration as 
specified by ispray_mod. It is primarily a controlling 
routine for spray computations. 

This routine has the following functions: 

(1) It initializes the source terms to zero. 

(2) Checks to see if new particles need to be introduced. 

(3) Advances liquid phase equations over an allowable or 
pre-specified time step, dtml, with calls to the following 
routines: 

intplal - Interpolates the gas phase properties at the particle 
location. 

chaslv - Advances liquid phase equations. 

intpla - Identifies computational cells and PEs associated 
with spray particles. 

sprips - Evaluates the liquid phase source term contributions 
of the CFD and PDF equations. 

(4) (a) For unsteady spray computations (steady _spray_model=F), 
this routine is called once at the beginning of every 

global time-step, dtgl. This is accomplished by continuing 
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Table 9. Description of LSPRAY-V Fortran subroutines & functions (continued) 


Subroutine 

Purpose of the Subroutines 

with steps (2) and (3) until the computations are completed 
over one single global time step, dtgl. (b) For steady spray 
spray computations (steady .spray _model=T), this routine is 
only called once to solve entire lifetime of all particles 
injected. It repeats step (3) until there are no more spray 
particles left to be integrated (note: The particles are taken 
of computation when they reach either a certain size of 
negligible proportion or exit out of the computational domain). 

(5) Returns control over to other solvers, e.g. CFD and EUPDF, 
and supply them with the source terms averaged over dtgl. 

dropdis(rhol, 
flowdum,sr, 
fid, Sind, nofg) 

This routine provides initial droplet distribution from 
the following correlation (also, appears as Eq. 91): 

dn/n = a{{D/D 32 ) al P)exp{-b((D/D 32 ) bet ))dD/D 32 
where a, b, alp, and bet are constants. 

drop_ic(n_i, 

nmih,nmis, 

nmip) 

It is called by the spray_main_injection routine in conjunction 
with spray _table(n J))=F at the beginning of every injection 
time step, dtil, for introducing a new group of spray particles. 

It applies for the introduction of new particles in conjunction with 
Options (2) or (3) of Sec. 9. In Opt. (2) the initial conditions 
are specified from the known correlations and in Opt. (3) from 
the available primary atomization models. 

find_transport_ds( 

ij le , chem _mo del , 

numb er_of .species, hp_eos, 

temp.gas , pres .gas ,y .gas , 

rhogm,visgm,congm,difgm) 

It computes the following properties of a gas mixture 
at the droplet interface by making use of the one-third 
rule of Eq. (46) : molecular viscosity, gas density, 
thermal conductivity, and diffusion coefficient. 

find jxyz face 

(i,xcfac,ycfac, 

zcfac) 

This routine computes x, y, and z locations of all the face 
centers of an element, i. This information is used in the 
particle search algorithm. 

get _liq_t at .properties 
(tmp,pmp,j) 

It computes the following variable properties of a liquid 
mixture: density, specific heat, molecular viscosity, 
gas density, thermal conductivity, and diffusion coefficient. 
This routine also computes the surface tension of a 
multi-component liquid fuel. 

get .liq.t at .properties 1 

(tmp,pmp,boil_tem, 

cp_sph,hvap_sph,icuin) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: boiling temperature, specific 
heat at constant pressure, surface tension, heat of vaporiza- 
tion, molecular viscosity, & liquid density and volume. 

get_liq_tat_properties_sp 
(tmpl ,tmp2,cplt 1 , 
cplt2,j) 

This subroutine is used in conjunction with the superheat 
vaporization model and computes the following variable 
properties of the liquid fuel: specific heat at droplet 
temperature, & specific heat at boiling temperature. 
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Table 9. Description of LSPRAY-V Fortran subroutines & functions (continued) 


Subroutine 

Purpose of the Subroutines 

init_psi2_dsd 

This routine provides initial droplet distribution from 
a clipped probability distribution function as given by 
Eq. 92. 

intpla 

This routine performs three main functions: 

(1) Particle Tracking - It identifies the computational cell 
in which a particle is located. In parallel computing, 

it also means identifying the corresponding processor of the 
computational domain in which a particle is located. 

(2) It implements appropriate boundary conditions. 

(3) It reassigns the particles between different PEs based on 
the parallel strategy employed. 

intplal 

This routine interpolates the gas-phase properties at the 
particle location. In the present case, a simple first-order 
interpolation is employed. 

niimd spray 

For computational elements where the neighboring cells are 
partitioned between different processors, this routine initializes 
the arrays, ipr_fr_id() and ile_fr_id(), by storing the information 
on the corresponding processor and element ID numbers of 
all the neighboring cells. This information is needed in order to 
track the particle movement between neighboring cells. 

mimd_spray_recv 

(i_recvfrom) 

This subroutine is called by mimd_spray in order to gather 
some relevant information from the neighboring processors. 

mimd_spray_send 

(Lsendto) 

This subroutine is called by mimd_spray in order to send 
some relevant information to the neighboring processors. 

par_loc(xparz, 

yparz,zparz, 

ipare,iparp) 

Given the x,y,z coordinates of a particle location, the 
algorithm identifies the corresponding computational cell 
in which a particle is located. It also identifies the 
corresponding processor of the computational grid in which 
a particle is located. This information is used to locate 
newly injected droplet groups. 

prnspr 

It provides written output of the spray data in a standard format. 

spray _int_rerun 

This routine has two functions: (1) It initializes the spray 
solver based on the data read from various input files. 

(2) If it is a rerun, it restarts the computations from the 
stored data of a previous computation. 
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Table 9. Description of LSPRAY-V Fortran subroutines & functions (continued) 


Subroutine 

Purpose of the Subroutines 

spray _main_injection() 

This routine is called to inject a new group of spray 
particles for each injector at the beginning of every 
injection time step, dtil. 

spray _output 

This routine is called to create appropriate restart files 
for the spray computations and also to provide for debugging 
purposes written output of some useful gas and spray data in 
a standard format. 

spray _plot_output 

It provides written output of some useful spray data in 
a standard format. It can be used in conjunction with 
steady _spray .model = F. 

spray cread-parameters 
(ipid,ncells)) 

This routine is called once in the beginning to read and 
initialize some of controlling parameters associated with 
the spray computations. 

spray .read-transport 
( (ipid, fid, filename) 

This routine is called once to read liquid thermo and 
tranport data: various constants used in the evaluation 
the polynomials in temperature for constant pressure 
specific heat, thermal conductivity, and viscosity. 

sprips 

This routine is called to provide the liquid-phase source terms of 
Eqs. (1)~(4) for use in both CFD and Monte Carlo PDF solvers: 

smlc(i) = liquid-phase contribution of Eq. (1) of Section 
2. 

smlmx(i), smlmy(i), smlmz(i) = liquid-phase contribution 
of Eq. (3) of Section 2. 

smle(i) = liquid-phase contribution of Eq. (4) of Section 2. 

sy(il,iu,bb, 

dd,aa,cc) 

It is a tri-diagonal matrix solver. It is used in the solution 
of both Eqs. (27) & (37). 

uvw.par (swlr_angle(nx) , 
angle .work , nx , ny, nn , 
v_inj ,nmip,t-rotation, 
cone,n_cone_rays, 
cone_rotation,uloc, 
vloc,wloc,n_l) 

It computes the initial particle injection velocity for different 
geometric configurations of spray as shown in Figs. 2 to 4. 
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Table 10. Description of PENG_ROB Fortran subroutines & functions 


Function 

Purpose of the Function 

peng_rob .density 
(chem_model,ycomp, 
tij, pressure, 
number _of_species) 

It provides a solution for the density of a multi-component 
mixture based on the Peng-Robinson equation of state. 

peng_rob .pressure 
(chem_model,ycomp, 
tij , density, 
number _of_species) 

It provides a solution for the pressure of a multi-component 
mixture based on the Peng-Robinson equation of state. 

Subroutine 

Purpose of the Subroutines 

hp.diff.table 

(ct,pt,dabcvi) 

This subroutine provides a value for P)/(D j 4 b P) + 

based on the Takahashi high-pressure correlation for the 
calculation of the diffusion coefficients in a gas mixture. 

It is based on the interpolation from the data derived from 
a graphical illustration of Reid et al [16], “The Properties 
Gases and Liquids.” The retrieved data is stored in 
hp_diff_table.dat as a table of f(tpr_diff,ppr_diff) values. 

hp .diffusion 
(tij ,pij,its, 
nwarn,diin_hp) 

It is called to evaluate the gas-phase mixture-averaged 
diffusion coefficients valid at high pressures. 

hp.transport 
(calc_name,its, 
tij TO ,nwarn, 
avisi_hp,condi_hp) 

It is called to evaluate the gas-phase transport properties 
valid at high pressures: viscosity and conductivity. 

peng_rob_all_rhos( 
chem_model,ycomp, 
tij, pressure, 
number _of_species, 
rhol, rho2, rho3, 
vvl, vv2, vv3, 
zl, z2, z3, 
aprgm,bprgm) 

It provides the roots for the cubic equation, Eq. (6), involving 
the compressibility factor, Z. One of the three roots of Eq. (6) 
yields Z„ and the other Z i. Jeff Moder made significant 
improvements to this routine. 

peng_rob_gen 

(xsp,ppr,tpr, 

its,rho0) 

It takes the following as the input variables: xsp() - the mass 
fraction of the species, ppr - pressure, tpr - temperature, and 
its - the number of species. And it returns as the output, rhoO - 
density based on the solution of the Peng-Robinson EOS. 

peng_rob_zcrit_rhor 

(ppr,tpr,its) 

It provides a solution for Z i C and pi r for ith species based 
on the solution of the Peng-Robinson EOS for a given value 
of P i r and Tj r . It also requires the specification of the 
total number of species, its. 
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Fig. 11a An illustration of the parallelization strategy 
employed in the gas flow computations. 
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Fig. lib An illustration of the parallelization strategy 
employed in the spray computations. 
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Table 11. CPU time (sec) per cycle versus number of PEs. 


Number of processors 

Solver 

Characteristic 

2 

5 

10 

CFD 

5 steps/cycle 

2.50 

1.25 

0.75 

EUPDF 

1 step/cycle 

6.5 

2.9 

1.9 

LSPRAY 

100 steps/cycle 

1.70 

0.64 

0.53 


Table 12. CPU Time (Sec) Per Cycle Versus Number of PEs. 


Number of Processors 

Solver 

Characteristic 

1 

2 

4 

8 

16 

Spray 

100 steps/cycle 

6.83 

5.29 

2.94 

1.64 

0.87 

Max. Spray Particles in a PE 

2695 

2097 

1165 

623 

312 

Min. Spray Particles in a PE 

2695 

598 

118 

14 

0 


• Both the EUPDF and CFD modules are well 
suited for parallel implementation. For the gas- 
phase computations, the domain of computation 
is simply divided into n-Parts of nearly equal 
size and each part is solved by a different pro- 
cessor. Fig. 11a illustrates a simple example 
of the domain decomposition strategy adopted 
for the gas-phase computations where the total 
domain is simply divided equally amongst the 
available computer processing elements (PEs). 
In this case, we assumed the number of available 
PEs to be equal to four. 

• But the spray computations are more difficult to 
parallelize for the reasons summarized below: 

(1) Non-uniform nature of spray distribution: 
Most of the particles are usually confined to a 
small region where the atomizer is located. 

(2) Dynamic nature of Lagrangian particles: 
Particles keep moving between different sub- 
divided domains of an Eulerian grid (grid used 
in the CFD computation) which are assigned to 
different processors. While some new particles 
are introduced at the time of fuel injection, some 
others are taken out of computation. 

Conceptually, there are several ways to paral- 


lelize the spray computations, we, however, devel- 
oped and tested two different domain decomposition 
strategies [10-14], 

• Strategy I: 

The Lagrangian particles were assigned fairly 
uniform amongst the available processors but the 
calculations associated with the particle track- 
ing, the interpolation of the gas-phase proper- 
ties, and the source-term evaluation were com- 
puted on the processor of the computational grid 
in which a particle is located. 

This strategy leads to an uniform loading during 
integration but leads to excessive message pass- 
ing. 

• Strategy II: 

The Lagrangian particles were assigned to the 
processor of the computational grid where the 
particle is located. Fig. lib illustrates a sim- 
ple example of the domain decomposition strat- 
egy adopted for the liquid-phase computations 
where the corresponding gas-flow computational 
domain is divided into equal parts between the 
four available PEs. In this strategy, the La- 
grangian particles are assigned to the processor 


N AS A/CR— 20 15-218918 


49 


of the computational grid where a particle is lo- 
cated. 

This strategy lead to a non-uniform loading dur- 
ing integration but leads to less message passing. 

Our experience has shown that Strategy II 
seems to work well on different computer platforms: 
both massively parallel computers as well as hetero- 
geneous cluster of workstations. So in the present ver- 
sion of the code, we have opted to implement Strategy 
II over Strategy I. 

19 DETAILS OF PARALLEL 
PERFORMANCE 

The details of the combined parallel perfor- 
mance of the CFD, EUPDF, and LSPRAY codes in- 
volving several different cases can be found in Refs. 
[10-15]. Here, we only summarize briefly the parallel- 
performance results for two different cases. One is 
a 3D test case and more details on this case can be 
found in the reference [13]. For this case, the calcu- 
lations were performed on a computational grid com- 
prising of 8430 tetrahedral elements and 100 Monte 
Carlo PDF particles per cell. The computations were 
performed on one of the NASA Ames Research Cen- 
ter’s parallel computer platforms called Turing which 
is a SGI Origin work-station with 24 PEs (Proces- 
sor Elements). Table 11 summarizes the CPU times 
per cycle taken by the EUPDF, LSPRAY, and CFD 
solvers vs the number of PEs. Both the CFD and 
PDF solvers show good parallel performance with an 
increase in the number of processors but for the spray 
solver it shows reasonable parallel performance. 

Next, we would like to summarize the results 
from [13, 38] showing only the results of spray compu- 
tations. The results are summarized in Table 12. The 
computations were performed on Turing at NASA 
ARC (Ames Research Center) - it is a SGI Origin 
work-station with a maximum of 24 PEs (Processor 
Elements). The computations made use of an un- 
structured grid with a mesh size of 3600 tetrahedral 
elements. And it made use of about 2,695 Lagrangian 
particles for the spray computations and one hundred 
Monte Carlo particles per element for the PDF com- 
putations. In a given cycle of one global time step, 
dtgi , the spray equations were advanced over one hun- 
dred time steps as given by dtgi/dt m i = 100. The 
first row of Table 12 summarizes the CPU times per 
PE per cycle taken by the spray solver vs the num- 
ber of PEs. As expected, the CPU time goes down 
with an increase in the number of processors. In an 


ideal case, one would expect an inverse reduction in 
epu time with an increase in the number of proces- 
sors. Here, we don’t get such an ideal performance 
because of the resulting non-uniform distribution of 
spray particles, between various participating proces- 
sors, from the implementation of Stratergy II. To get 
an idea of the spray particle distribution, we have tab- 
ulated the maximum and minimum number of parti- 
cles found between various processors. When we go 
from 1 to 2 PEs, 2097 particles are assigned to one 
and the rest to the second. With four they are dis- 
tributed between 1165 and 118, with eight between 
623 and 14, and with sixteen from 312 to 0. The 
results clearly show that the reduction in the CPU 
time varies almost linearly with the reduction in the 
number of maximum particles. 

20 SUMMARY OF SOME VALIDATION 

CASES INVOLVING BOTH REACTING 
AND NON-REACTING SPRAY 
COMPUTATIONS 

A total of the following nine cases were vali- 
dated: 

1. A reacting methanol spray with no-swirl, 

2. A non-reacting methanol spray with no-swirl, 

3. A confined swirl-stabilized n-heptane reacting 
spray, 

4. An unconfined swirl-stabilized n-heptane react- 
ing spray, 

5. A confined swirl-stabilized kerosene reacting 
spray, 

6. A flashing jet generated by the sudden release of 
pressurized R134A from cylindrical nozzle, 

7. A liquid jet atomizing in a subsonic cross flow, 

8. A Parker-Hannihn pressure swirl atomizer, & 

9. A single-element LDI (Lean Direct Injector) 
combustor experiment. 

The experimental data for the first two cases 
was provided by McDonell & Samuelsen from the 
University of California at Irvine [38]. Both the cases 
are without swirl; one is a reacting case and the other 
is non-reacting. The data for the third and fourth 
cases was provided by Bulzan from the NASA Glenn 
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21 CONCLUDING REMARKS 


Research Center [39-40]. Both the cases are swirl- 
stabilized reacting cases, one is an unconfined flame 
and the other is confined. The data for the fifth case 
was provided by El Banhawy & Whitelaw from Im- 
perial College [4]. It is a confined swirl-stabilized 
kerosene spray flame. Further details of the last four 
cases can be found in [78]. These cases were chosen 
because of their importance in some aerospace appli- 
cations. 

Here, we provide a brief summary of our val- 
idation effort but a detailed presentation of the re- 
sults and discussion can be found elsewhere in the 
papers [10-13, 78]. The validation is based on some 
3D and axisymmetric calculations involving both re- 
acting and non-reacting sprays. Some of them were 
performed on unstructured grids and the others on 
structured grids. 

The comparisons involved both gas and drop ve- 
locities, drop size distributions, drop spreading rates, 
and gas temperatures. In general, the predicted re- 
sults provide reasonable agreement for both mean 
droplet sizes (Z? 32 ) and average droplet velocities but 
mostly underestimate the droplets sizes in the inner 
radial region of a cylindrical jet. 

The comparisons also involved the results ob- 
tained from the use of the Monte Carlo PDF method 
as well as those obtained from a conventional CFD 
solution without the Monte Carlo PDF method. For 
the first case of McDonell & Samuelsen’s reacting 
spray flame, the detailed comparisons clearly high- 
lighted the importance of chemistry /turbulence in- 
teractions in the modeling of reacting sprays [13]. 
The results from the PDF and non-PDF methods 
were found to be markedly different with the PDF 
solution providing a better approximation to the re- 
ported experimental data. The PDF solution showed 
that most of the combustion occurred in a predom- 
inantly diffusion-type of flame environment and the 
rest occurring in a predominantly premixed-type of 
flame environment. However, the non-PDF predic- 
tions showed incorrectly that most of the combustion 
occurred in a predominantly vaporization-controlled 
regime. The Monte Carlo temperature distribution 
showed that the functional form of the PDF for the 
temperature fluctuations varied substantially from 
point to point. The results brought to the fore some 
of the deficiencies associated with the use of assumed- 
shape PDF methods in spray computations. 


• This manual provides a complete description of 
LSPRAY-V - a Lagrangian spray solver devel- 
oped for application with parallel computing and 
unstructured grids. 

• It facilitates the calculation of the multi- 
component liquid sprays with variable properties 
valid over a wide range of pressure conditions en- 
countered in gas-turbine engines. 

• It provides the user with a basic understanding 
of the the spray formulation and the LSPRAY-V 
code structure, and complete details on how to 
couple the spray code to any other flow code. 

• The basic structure adopted for the grid repre- 
sentation and parallelization for the gas side of 
the flow computations follows the guidelines es- 
tablished for NCC. 

• Also, we have extended the joint scalar Monte 
Carlo PDF method to two-phase flows and, 
thereby, demonstrating the importance of chem- 
istry/turbulence interactions in the modeling of 
reacting sprays. 

• Based on the validation studies involving several 
confined and unconfined spray flames, the results 
were found to be encouraging in terms of their 
ability to capture the overall structure of a spray 
flame. 

• The source code of LSPRAY-V will be available 
with OpenNCC as a complete package. 
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APPENDIX 

DETAILS OF THE PRIMARY 
ATOMIZATION MODELS 

1. Sheet Breakup Primary Atomization 
Model (Likely Application: Pressure Swirl 
Atomizer) 

Here, we summarize the details of the sheet 
breakup model taken from Schmidt et al [48]. The 
growth of an infinitesimal disturbance as given by 

77 = r] 0 exp(ikx + cut) ( 1 ) 

was analyzed based on a linear stability analysis of a 
two-dimensional, viscous, incompressible liquid sheet 
of thickness 2h which moves through an inviscid, in- 
compressible gas medium. This was analyzed in a co- 
ordinate system moving with the sheet with a relative 
velocity of U where r] 0 is the initial wave amplitude, 
k (= 2tt/X) is the wave number, and ui = uj r + iuii 
is the complex growth rate. The most unstable dis- 
turbance responsible for the sheet breakup is denoted 
by H. 

Based on the linearized liquid and gas continu- 
ity and momentum equations subject to the linearized 
boundary conditions at the gas and liquid interfaces, 
a sinuous mode dispersion relation was obtained by 

[49], 


cu 2 [tanh(kh) + Q] + u)[Ai'ik 2 tanh(kh) + 2iQkU]-\- 
Av 2 k A tanh(kh ) — 4 i/fk 3 l tanh(lk) — QU 2 k 2 + — 0(2) 


where Q = p g /pi, l 2 = k 2 + w/vi, and U is the 
relative velocity between the liquid and gas. Inviscid 
analysis also indicates that for low gas Weber num- 
ber (We = p g hU 2 /a) flows, long waves tend to grow 
leading to liquid sheet breakup but for higher Weber 
numbers, short waves produce a maximum growth 
rate followed by breakup. The critical Weber number 
that leads to the transition from the long wavelength 
regime to the short wavelength regime was shown to 
be We g = 27/16. For most modern fuel injection 
systems, the film Weber number is well above this 
critical limit. The growth rate for the sinuous mode, 
ui r , based on an order of magnitude analysis of the 
dispersion relation yields, 

2 ui k 2 1 anh(kh) 

CJ r - 1 “ 

tanh(kh) + Q 


4i k^tanh 2 (kh) — Q 2 U 2 k 2 — [tanh(kh) + Q](—QU 2 k 2 + crk 3 /pi ) 


tanh(kh) + Q 


( 3 ) 


For short waves in the limit of Q << 1 for high-speed 
sheets, it yields 


U) r 


-2 vik 2 + a Uuf fc 4 + QU 2 k 2 


ak 3 

Pi 


( 4 ) 


Following Dombrowski & Johns [51], the sheet 
disintegration leads to the formation of ligaments 
once the unstable waves reach a critical amplitude 
and Eq. (4) shows that the growth rate of short 
waves is independent of the sheet thickness. The cor- 
responding breakup time r and the breakup length L 
are given by: 


? 7 b = t] 0 exp(Hr) => r = In — (5) 

if rj 0 

L = Vt = ^ ln( ( 6 ) 

where H is the maximum growth rate as determined 
by Eq. (4), the term ln( — ) has an assigned value of 
12 as suggested by Dombrowski & Hooper [50], and 
V is the absolute velocity of the liquid. 

The initial diameter of the ligaments is derived 
from a mass balance relationship. For long waves, it 
is assumed that the ligaments are formed from tears 
in the sheet once per wavelength and the resulting 
diameter is given by, 


dL = 



( 7 ) 


where K s b is the wave number corresponding to the 
maximum growth rate fl as obtained from Eq. (3) 
and the film thickness, h, is calculated from the 
breakup length, L, the radial distance from the cen- 
terline to the mid-line of the liquid sheet at the at- 
omizer exit, r Q , and the spray angle, 0 , as follows: 
h = 2, Pl v{r 0 +L^n(e/2)) - For short waves, the lig- 

ament diameter is independent of the liquid sheet 
thickness and is assumed to be proportional to the 
wave length associated with the maximum growth 
rate f 1 as follows: dL = 2 )/// where the ligament 

constant, Cl, is equal to 0.5. 

For both long and short waves, Dombrowski & 
Johns [51] developed a linear stability analysis for 
the further breakup from ligaments to droplets based 


N AS A/CR— 20 15-218918 


52 


on the Weber’s analysis on capillary instability. The 
analysis shows that the breakup occurs when the 
amplitude of the unstable waves nears the radius of 
the ligament. And the corresponding most unstable 
wavenumber, AAb, is given by: 


K^bdL 


3 pi 


2 yj piadi) 


(8) 


This analysis thus yields the most probable droplet 
size based on a simple mass balance calculation where 
d? b = 37 xd 2 L /K Lb . 


1.1 Application to pressure swirl atomization 


For the pressure swirl atomizer, the initial in- 
jector exit velocity and liquid sheet thickness are cal- 
culated following the approach of Schmidt et al [48]. 
It assumes a uniform velocity profile for the initial 
liquid velocity, V, as given by, 


V = max{ 0.7, 


4m 

7 idgPicosO 



( 9 ) 


where m and 6 are the measured mass flow rate and 
spray angle, respectively, d 0 is the injector exit di- 
ameter, and A p is the pressure drop in the injector. 
Once V is known, the corresponding axial component 
of the sheet velocity is calculated via u = V cosd. 
And the initial sheet thickness h a is calculated from 
the conservation of mass: 


m = Trpiuh 0 (d 0 - h 0 ) (10) 

At the point of primary breakup, the actual drop size 
is chosen from a Rosin-Rammler distribution with the 
mean size as given by db of the sheet breakup model. 
Further movement of the droplets is tracked by mak- 
ing use of a Lagrangian formulation. 

2. Blob Jet Primary Atomization Model 
(Likely Application: Single-Orifice Nozzles) 

Here we summarize the details of the blob 
jet primary atomization model taken from Reitz & 
Bracco [41] and Reitz [52,45]. It applies for a cylin- 
drical liquid jet issuing into an incompressible gas. 
The following dispersion relation was obtained based 
on the stability analysis of a cylindrical liquid surface 
subjected to linear perturbations: 


2 _ r, ,2 r J i( fca ) 2kl h(ka) I((la) 1 
w w i /o ( fca ) k 2 + l 2 I 0 (ka)I 0 (la) i 


ak n 

pia 2 \ n P + k 2i I 0 (ka) 

pi l z + k z lo\ka) ARfca) 


( 11 ) 


where io, A, and Ao, A'i are the modified Bessel 
functions of the first and the second kinds. 

Reitz [52,45] generated curve-fits of numerical 
solutions to Eq. (11) for the maximum growth rate 
(w = fi) and the corresponding wavelength (A = A): 


A 

a 




(1 + 0.45Z°' 5 )(1 + 0.4T 0,7 ) 
(1 + 0.87Wei'67)0.6 

0.34 + 0.38 Wef 5 
(1 + Z)( 1 + 1.4T 06 ) 


( 12 ) 

(13) 


where Z = , T = ZWe° g 5 , Wei = We g = 

Pg ^ a , and Rei = . A core region is predicted with 

the blob injection method because there is a region of 
large discrete liquid parcels near the nozzle. Based on 
the jet stability theory, new drops are formed from a 
parent drop or blob. It is assumed that small droplets 
(with radius, r) are formed from the parent drops 
(with radius, a) with drop size proportional to the 
wavelength of the fastest-growing or most-unstable 
wave, 

r = B 0 A (if B 0 A < a) or 

r = mm{(37ra 2 t//2H) 0 - 33 , (3a 2 A/4) 0 - 33 } 

(if B 0 A > a, one time only ) (14) 

where B a = 0.61 according to Reitz [52,45]. In the 
above, it is assumed for the (B a A < a) condition 
that small droplets are formed with the dropsize pro- 
portional to the wavelength of the fastest growing 
mode and the second (B a A > a) condition applies to 
drops larger than the jet and it assumes that the jet 
disturbance has frequency f2/ 2ir ( a drop is formed 
each wave period) or that the drop size is determined 
from the volume of liquid contained under one sur- 
face wave. And the rate of change of droplet radius 
due to breakup is given by 


da/dt = — (a — r)/r (r < a) (15) 

where r is the breakup time, r = 3.726Ria/AH, and 
the value for the breakup time constant, B i, depends 
on the injector characteristics and its value ranges 
between 1.732 to 40. And aft = t a ) = a Q is the 
initial drop radius at time, t 0 . 
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After the breakup, a new parcel containing 
product drops of size, r, is created and added to the 
computations. This was done as long as the mass of 
the liquid removed from the parent (p;47r(a 3 — a 3 )/3) 
reached or exceeded 3% of the average injected parcel 
mass and if the number of product drops is greater 
than or equal to the number of parent drops [45]. 
While waiting for sufficient product drops to accu- 
mulate, the parent drop number was adjusted so that 
Na 3 = N 0 a 3 but the parent drop number, N 0 , was 
then restored following the creation of the new prod- 
uct parcel [45]. 

In the case of (B a A > a), the parent parcel 
was replaced by a new parcel containing drops with 
size given by Eq. 14 after a time equal to r (with 
N = N 0 a 3 /r 3 ) [45]. This breakup procedure was 
allowed only once for each injected parcel [45]. 

Validation of the model for a single hole orifice 
in a typical diesel engine was demonstrated by Reitz 
and Diwaker [46] and Reitz [52,45]. 


Here we summarize the details of the BLS at- 
omization model taken from Khosla and Crocker [47] . 
In this model both surface shear breakup and col- 
umn breakup modes are included. Before column 
breakup, fragments may be formed due to boundary 
layer stripping depending on the local Weber num- 
ber and q (= piuf / p g u 2 g ). When the jet reaches the 
column breakup time, the entire jet breaks into frag- 
ments. It also allows for further breakup of the frag- 
ments based on a modified boundary layer stripping. 
And it is followed by a final breakup step based on the 
Rayleigh-Taylor secondary droplet breakup model. 

Fragments are stripped from the liquid column 
if the We g satisfies the following criteria: 

We g > 50 Re] 12 g _1/0 - 81 and 

We g > 15 (18) 

where 


3. Air Blast Primary Atomization Model 
(Likely Application: Air Blast Atomizers) 

The air blast atomization model is essentially 
based on the idea of pressure-swirl atomization model 
(Section 1.1) as the primary atomization of an air 
blast injector is based on the aerodynamic analysis 
involving the Kelvin-Helmholtz instability of a liq- 
uid jet in an incompressible gas. It, however, differs 
from the pressure-swirl atomization model in the de- 
termination of the initial sheet velocity and thickness, 
which are given by 


V s heet = aVi + (l-a)V g (16) 


S 


r [l 



(17) 


where a has a value of 0.12 to 1 depending on the 
fuel filmer characteristics, r is the radius of the fuel 
filmer, and to is the fuel flow rate. The continuous 
annular sheet is represented by a finite number of 
point injectors located randomly along the circular 
ring of the liquid sheet. 


4. Modified BLS Primary Atomization 
Model (Likely Application: Liquid Jet in a 
Cross Flow) 


Re g — PgdjUg/ p g 
we g = u 2 g djp g /ai 

where Re g and We g are based on the gas velocity 
component normal to the liquid jet direction instead 
of the relative velocity. If column stripping does oc- 
cur, the amount of mass removed from the column is 
given by 


M s hed — 


3 

4 


ndpi^-u re iAa 

t 



where, 


rPy I V 3 r/fff I !/ 3 
P; Pi 

r 8 /q -1 1/2 

3 Au re ipi 



(19) 


(20) 

(21) 


(22) 


where tb is the liquid column drop lifetime. The ad- 
dition of the factor t* /tb causes the shedding rate to 
increase essentially linearly with distance away from 
the injection location. This accounts for the lack of 
shedding close to the injection location and the sub- 
sequent buildup of shedding over the life of the liquid 
column. The shed drop SMD (Sauter Mean Diame- 
ter) is given by 
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(23) 


the same as Eqs. (24) and (25). The lateral velocity 
component is given by 


SMD = 3.1— d 1 / 2 [—1 1/A [—tH — 1 1/2 

t* L p g 1 L U re lpi 

The t b /t* factor is included with the effect of 
producing smaller drops near the injection location. 
However, tb/t* is limited to a minimum value of 2.5 
which is never exceeded for some cases. The amount 
of mass shed is tracked and 10 new parcels are cre- 
ated when the cumulative shed mass after a time step 
exceeds 1% of the mass of the parcel. Each parcel is 
allocated an equal amount of the shed mass and the 
size for each new parcel is selected randomly from a 
uniform distribution between 0.4 and 1.6 times the 
mass mean diameter, MMD. And the drop velocities 
were given by 

Ud = Up + 0.3(RND — 0.2) (u g — u p ) (24) 


Vd = v p + +0.3(RND — 0.2)(v g — v p ) (25) 


Wd = Wp + 0.25 (RND — 0.5 )(u re i — w p ) (26) 


Wd = Wp + 0.1 (RND — 0.5) [u re i — w p ) (28) 

Fragments are further broken into small drops 
according to a modified version of the boundary layer 
stripping model based on the following criteria, 

We g > \/~Re and 

We g > 15 (29) 

where 

we g = u 2 rel d d p g /(Ji 

R&g — PgddU re l [ p g 

The criteria are generally the same as for column 
stripping except that the dependence on q is not 
needed. Also note that We and Re are now deter- 
mined using the relative velocity instead of the cross 
flow velocity. Again, the mass shed from a fragment 
in a time step and the SMD are given by 


where RND is a random number. Column stripping 
occurs, assuming the above criteria are met, until the 
column breakup time is exceeded. Parcels created 
through the column stripping mechanism are con- 
sidered to be fragments which may undergo further 
breakup as discussed below. First, though, the col- 
umn breakup mechanism, which also produces frag- 
ments, is described. The liquid jet column breakup 
time is given by 


t b = A b We 0 62 t* (27) 

The constant, A b , has a value of 25. Since the present 
model includes a fragment breakup step after the col- 
umn breakup, the We dependence for the breakup 
onset was retained. 

After the column breakup time is reached, the 
column is broken into 18 new parcels with MMD = 
0.45dj. The new parcels are also designated as frag- 
ments. The size of the fragments still tend to be 
large, so the ultimate drop size from the primary 
breakup process is mostly determined from the frag- 
ment breakup process. The size distribution is the 
same as described above for the column stripping. 
The cross flow and normal velocity components are 


M s hed 


1 . 27 T dpiu re iAa 



(30) 


SMD = 3.6d 1/2 [—1 1/4 ^ — 1 1/2 (31) 

where A and a are given by Eqs. (20) and (21), re- 
spectively. The new droplet velocities are given by 
Eqs. (24), (25), and (28). The broken fragments pro- 
duce 3 new parcels with size distribution the same as 
described above when the shed mass from the frag- 
ment exceeds 20% of the fragment mass. A fragment 
can continue to breakup until it no longer meets the 
criteria of Eq. (29) or until its size is lower than the 
newly created drops. Once the fragment breakup pro- 
cess is complete, drops may breakup further based on 
a Rayleigh-Taylor secondary breakup method. 

Khosla and Crocker [47] applied the model to 
predict the properties of Jet A-l kerosene fuel in- 
jected into a cross-flowing air stream. 
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DETAILS OF THE SECONDARY 
DROPLET BREAKUP MODELS 

Recent studies by Reitz et al [55, 52] have ex- 
amined the breakup of single droplets moving in 
a transverse, high-velocity air jet. The high-speed 
photography provided new insights into the details 
of the breakup mechanism of a single drop. The 
droplet breakup regimes are classified as bag, strip- 
ping (shear) and catastrophic (surface wave) based 
on an increasing size of Weber number. In the bag 
breakup mode (at low Weber number), the drop is 
flattened by the aerodynamic pressure, then turned 
inside out, forming a thin hollow bag which is tied to- 
gether with a circular belt-like structure on the wind- 
ward side. The bag eventually bursts into smaller 
liquid fragments, whereas the belt decays into larger 
ligaments and droplets. In the stripping regime thin 
sheets or ligaments of fluid are continuously shed from 
the periphery of the distorting parent drop as a con- 
sequence of a K-H instability, causing these sheets to 
disintegrate into tiny droplets. This process always 
leaves a coherent residual parent drop. The catas- 
trophic drop breakup takes place in two stages lead- 
ing to a collection of larger and tiny product droplets: 
Large amplitude long-wavelength waves caused by 
drop deceleration induce a R-T instability on the 
flattened drop which leads to a breakup into large 
product droplets, while at the same time short sur- 
face waves induce a K-H instability on the windward 
side of the parent drop resulting in a collection of 
much smaller product droplets. In diesel or other 
high-pressure gas-turbine sprays the droplets span a 
wide range of velocities and hence Weber numbers, 
and thus it is expected that all three droplet breakup 
mechanisms are relevant in the breakup modeling. 

In what follows, we provide some details of 
the secondary droplet breakup models contained in 
the CFDRC/UW atomization module: (1) Raylcigh- 
Taylor, (2) TAB, and (3) ETAB. These details are 
taken from [52-56, 42-44]. 

1. Rayleigh- Taylor Secondary Droplet 
Breakup Model 

Here we summarize the Rayleigh-Taylor sec- 
ondary breakup model developed by Patterson and 
Reitz [57]. It is based on the analysis developed 
by Taylor [53,57] that accounts for the disturbances 
caused by droplet deceleration. In the Rayleigh- 
Taylor breakup mechanism, the breakup wavelength, 
A, is given by 


A = 0.27ry // 3<r/|'Ud|(p/ - p g ) (32) 


where iid is the drop deceleration (= f Cd p 9 l d ■ c d is 
the aerodynamic drag coefficient, U is the drop rela- 
tive velocity, & d is the drop diameter). Furthermore, 
the breakup wavelength A is limited by 


A = max(0.8d, A) [1 + 0.2(RND — 0.5)] (33) 

and the breakup time, tb,RT, is given by 


tb,RT = Cfreq * 0.5y/a(pi + Pg)( 


\u d \{pi- Pg) 


1.5 


(34) 

However, the value assigned for the constant, Cf req , 
depends on the droplet classification - parent, prod- 
uct, or default. For more details, one can refer to 
Patterson and Reitz [21]. After the breakup, no new 
drop parcels are created and there is no change in ve- 
locities between the parent and product drops. How- 
ever, the drop number in a given parcel changes to 
R product as given by Rparent {^parent/ A) due to the 
change in the sizes between the parent and product 
drops. 


2. The TAB Secondary Droplet Breakup 
Model (Likely Application: 

Lenticular-Shaped Droplet Deformations) 

In an attempt to provide a description of the 
droplet and jet disintegration in the modeling of 
diesel sprays, O’Rourke & Amsden introduced the 
TAB model [42] . Here we summarize the TAB model 
taken from [43-44]. 

The TAB breakup model is based on Taylor’s 
analogy between an oscillating, distorting drop and a 
spring-mass system. A detailed analysis of this model 
together with a discussion of its numerical implemen- 
tation can be found in [42 & 58]. In this model, the 
drop motion is governed by a linear ordinary differ- 
ential equation for a forced, damped harmonic oscil- 
lator. The forcing term is given by the aerodynamic 
droplet-gas interaction, the damping is due to the 
liquid viscosity and the restoring force is supplied by 
the surface tension. The parameters and constants 
have been determined partly from theoretical consid- 
erations and partly from experimental observations. 

The TAB model describes the distortion of the 
drop by the deformation parameter, y = 2x/a , 
where x denotes the increase in the radius increase 
of the equator from its equilibrium position and a 
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is the drop radius. The equation for the distortion 
parameter y is given by 


V + 


5 fH 

Pia 2 


V + 



2p g \U\ 2 

3 pia 2 


(35) 


Assuming a constant relative drop-gas velocity, 
U (which is satisfied in the numerical solution process 
during a given time step), the solution to Eq. (35) is 
given by 


. , We 

«(t) = 


=i / r We-, 

e ‘d ( [y(o) — J cosut 




ut d 


sinuj 


l ) 


(36) 


Here we provide some details of the ETAB 
model taken from Tanner [43-44]. The ETAB model 
is based on the following modifications to the stan- 
dard TAB model: (1) the droplet disintegration is 
modeled via an exponential law which relates the 
mean product droplet size to the breakup time of the 
parent drop; and (2) an energy balance consideration 
between the parent and product droplets yields an 
expression for the normal velocity component of the 
product droplet. 

When the breakup condition of We > 
We cr it = 6 and y(t) = 1 is met, then the parent 
drop breaks up into a collection of product droplets, 
subject to a size distribution function which, in gen- 
eral, depends on the breakup mechanism. In the 
ETAB model, the rate of product droplet generation, 
dn(t)/dt, is given by 


where We = p g aU 2 /a and 


td 


2pia 2 
5 pi 


(37) 
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8a 1 

Pia 3 t 2 


(38) 


In this model, it is assumed that a necessary 
condition for drop breakup is reached when We > 
We cr it . And the value for the critical Weber num- 
ber is determined experimentally to be 6. For an 
inviscid liquid with initial deformation conditions 
2/(0) = y(0) = 0, the solution to Eq. (35) 

reduces to y(t) = We(l — cosw 0 f)/12, where 
uj 2 = [44]. Consequently, the breakup occurs 

when y(t) > 1. The drop size after breakup is de- 
termined by an energy balance equation between the 
parent and the product droplets which equates the 
surface energies of parent drops with the energies of 
product drops due to oscillation and distortion. Also, 
the product droplets are initially equipped with the 
additional deformation velocity x = ay/2, which 

acts normal to the path of the parent drop and is 
responsible for the formation of the spray angle. 

One major advantage of this model is that it is 
based on a simple linear equation and it can be used 
effectively to describe the lenticular-shaped droplet 
deformations as observed in the experiments of [56 & 
54], 


3. The ETAB Secondary Droplet Breakup 
Model (Likely Application: High-Pressure 
Diesel Engine) 


dJ f - **■*> 


(39) 


where n(t) = mo/fh{t) and mo is the mass of the 
parent drop and to the mean mass of the product 
droplet distribution. Utilizing the fact that dn/dt = 
— (mo/fh 2 )(drh/dt), leads to the breakup law which 
relates the product drop size to the breakup time as 
determined by the TAB model. 


— = -3 K b m (40) 

The breakup constant K b depends on the 
breakup regime and is given by parent drop prop- 
erties only. Bag breakup occurs if We = Wet and 
stripping breakup happens if We > Wet ■ And it is 
given by 


Kb = k\U> if We < We t or 

K b = k 2 bjVWe if We > We t (41) 

The values for We*, k\ and k 2 have been determined 
experimentally and has been set to fci « k 2 = 1/4.5 
and We t = 80. 

In this model, a uniform product droplet size 
distribution is assumed. It is also noted that the 
choice of uniform distribution is not expected to be 
realistic but may produce good approximations when 
averaged over many drop breakups, because parent 
drops of different sizes and Weber numbers will in 
general yield a wide range of duct droplet sizes. With 
this assumption, Eq. (40) becomes 

- = e~ Kbt (42) 

a 
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where a and r are the radii of the parent and product 
drops, respectively. 

After breakup of a parent drop, the initial de- 
formation parameters of the product droplets are set 
to y(0) = 2/(0) = 0. Also, the product droplets are 
initially supplied with a velocity component perpen- 
dicular to the path of the parent drop with a value 
vt = Ax, where A is a constant determined from the 
following energy balance consideration. The energy 
of the parent drop is the sum of the surface tension en- 
ergy and the droplet deformation energy. The second 
is computed as the product of the aerodynamic drag 
and the drop deformation at the stagnation point, 
estimated to be 5a/9. This leads to 


Eparent = 47rcra 2 + 5Trc d p g a 3 \ U\ 2 /18 (43) 

And the energy of the product droplets in the frame 
of reference of the parent drop is given by 


E P roduct = 4Trcra 3 /rsMR + A 2 tt pia 5 y 2 /6 (44) 

where the Sauter mean radius, vsmr , enters via the 
relation r 2 = a 3 /vsmr- From Eqs. (43) and (44) 
one obtains the relation 


A 2 = 3[1 - a/rsMR + 5c d We/72]u; 2 /y 2 (45) 

where w 2 = ■ 

Tanner [43-44] analysis yields an approximate 
value of 0.69 for A showing that only 70% of the par- 
ent drop deformation velocity goes into the normal 
velocity component of the product droplets, where as 
it is 100% in the standard TAB model. 

Also, the characteristic time, r (= A~), for 
breakup in Eq. (42) for an inviscid liquid {pi = 

0) is given by cniy^P' if We < We t or 
if We > Wet, where the suggested val- 
ues are for a\ = (v^fci)" 1 and «2 = (V8k 2 )~ 1 . 

The application of the ETAB model in the sim- 
ulation of a high pressure liquid jet breakup can be 
found in [43-44]. 
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